March 2021 xx xx 10.6339/21-JDS1010 December, 2020 February, 2021
11affiliationtext: Department of Applied Mathematics and Statistics, Colorado
School of Mines, Golden, CO, USA 22affiliationtext: Department of Applied
Physics and Applied Mathematics, Columbia University, New York City, NY, USA
33affiliationtext: Department of Engineering Sciences and Applied Mathematics,
Northwestern University, Evanston, IL, USA

# Investigating the Relationship Between Air Quality and COVID-19 Transmission

Laura Albrecht111All authors contributed equally to this work Paulina
Czarnecki Corresponding author. Email<EMAIL_ADDRESS>Bennet Sakelaris

###### Abstract

It is hypothesized that short-term exposure to air pollution may influence the
transmission of aerosolized pathogens such as COVID-19. We used data from 23
provinces in Italy to build a generalized additive model to investigate the
association between the effective reproductive number of the disease and air
quality while controlling for ambient environmental variables and changes in
human mobility. The model finds that there is a positive, nonlinear
relationship between the density of particulate matter in the air and COVID-19
transmission, which is in alignment with similar studies on other respiratory
illnesses.

###### keywords:

Air pollution; Effective reproduction number; Generalized additive model;
Mobility; Particulate matter.

##  1 Introduction

Poor air quality is a well-documented public health concern, causing
respiratory illnesses, negatively impacting cardiovascular health, and
decreasing life expectancy (Delfino et al., 2005; Krewski, 2009; Lelieveld et
al., 2015; Janssen et al., 2013). Especially harmful is fine particulate
matter (with diameter less than 2.5 $\mu$m, called PM2.5) suspended in the air
(Delfino et al., 2005; Laden et al., 2000; Schwartz et al., 1996). These
particles, which are mainly created by industrial combustion processes and
atmospheric reactions, can deeply penetrate the lungs and cause damage to the
respiratory system over time (Xing et al., 2016; Tucker, 2000). Furthermore,
it is hypothesized that viruses may attach themselves to these suspended
particles and infect hosts upon inhalation (Ciencewicki and Jaspers, 2007;
Liang et al., 2014; Contini and Costabile, 2020; Jiang et al., 2020). This
effect is particularly pronounced in cities with extreme levels of air
pollution, such as Beijing, where cases of influenza have been shown to
increase with the density of PM2.5 in the air (Ciencewicki and Jaspers, 2007;
Liang et al., 2014; Feng et al., 2016; Su et al., 2019; Jiang et al., 2020).
Since December 2019, when the novel coronavirus SARS-CoV-2 was first
identified in Wuhan, China, it has erupted into a global pandemic (Zhu et al.,
2020a; Jiang et al., 2020; Zhu et al., 2020b). Many risk factors and
comorbidities have been found, including obesity, pre-existing health
conditions such as diabetes or hypertension, and advanced age (Richardson et
al., 2020; Yang et al., 2020). An increase in environmental factors such as
air pollution, temperature, and humidity have also been found to be associated
with increased COVID-19 case counts early in the pandemic in China (Zhu et
al., 2020b; Jiang et al., 2020), however the effects of PM2.5 remain largely
unexplored as the pandemic continues to develop. Furthermore, recent events
such as raging wildfires in the western United States have increased
concentrations of PM2.5 in the air (O’Dell et al., 2019; Jaffe et al., 2008);
it would be useful to understand whether these changes can affect the spread
of COVID-19.
In this paper, we analyze the effective reproductive number ($R_{t}$) of the
disease, which is the average number of secondary infections caused by an
infectious individual on day $t$. Because of this time dependency, $R_{t}$
acts as a measurement of daily transmission intensity and as a result, is
frequently utilized for decision-making during a pandemic (Lipsitch et al.,
2011). After using the open-source software rt.live (Systrom et al., 2020) to
calculate $R_{t}$ in 23 provinces in Italy, we studied its association with
PM2.5 using a generalized additive models (GAM). In the GAM, we also
controlled for other environmental variables, including humidity, temperature,
and changes in mobility due to lockdown. We found that at concentrations of
PM2.5 greater than about 70 $\mu$g/m3, poor air quality is positively
correlated with increased $R_{t}$ of COVID-19. Furthermore, we performed a
23-fold cross-validation procedure in order to test the sensitivity of the
model on each province and test for selection bias. We found that our model
was robust against this test.

##  2 Methods

###  2.1 Study Area

In order to ensure that we had sufficient data to extract meaningful results,
we considered only Italian provinces that reported more than 2000 cumulative
cases of COVID-19 during the study period of February 24, 2020 to August 1,
2020. We then removed locations where environmental data was unavailable,
resulting in the final study area of the 23 provinces shown in Figure 1. These
provinces accounted for over 63% of total Italy’s COVID-19 cases over the time
frame.

Figure 1: A geospatial heat map with locations of the 23 provinces included in
the study, colored by the cumulative number of COVID-19 cases as of August 1,
2020.

###  2.2 Data

The environmental variables included in the study were daily median
temperature, relative humidity, and PM2.5. Temperature and relative humidity
were collected from the NOAA Integrated Surface Database using the worldmet
package in R (Carslaw, 2020), while PM2.5 data was taken from the World Air
Quality Index (https://aqicn.org/).

We also controlled for the effects of lockdown by incorporating data from
province-specific Google Community Mobility Reports
(https://google.com/covid19/mobility/), which quantify the daily change in
aggregate mobility compared to a pre-pandemic baseline. The dataset was
generated using anonymized data from Google Maps to measure daily visitors in
six categories of locations: residential areas, retail locations, transit
stations, grocery stores, parks, and workplaces. In order to gain perspective
on how mobility has changed due to the pandemic, these measures were divided
by a median value of visitors for that location for that day of the week,
which was measured prior the onset of the pandemic. For this study, we took
our mobility variable to be the change in mobility in residential areas;
however, we found similar results when we substituted data from any of the
other five categories of locations provided in the Google Community Mobility
Reports.

Daily tests counts and confirmed cases were collected directly from the
Italian Department of Civil Protection (https://github.com/pcm-dpc/COVID-19)
and were used to estimate the value of the effective reproductive number,
$R_{t}$.

$R_{t}$ was estimated using an open source model from rtlive2020 called
rt.live, which we adapted to be compatible with our data. Notably, the model
estimates total cases by scaling confirmed cases according to testing volume.
It also approximates a time series of cases by day of infection using an onset
delay distribution which was empirically driven by data from onset_delay while
assuming a five day incubation period as evidenced by work from incubation.
Finally, it uses the generation time derived by generation_time with mean 4.7
days and standard deviation 2.9 days to estimate $R_{t}$.

The time series of new COVID-19 cases, mobility, PM2.5, median temperature,
median humidity, and $R_{t}$ from February 24, 2020 to August 1st, 2020 are
shown in Figure 2. The peak in cases occurred between March 14th and April
18th in each province then cases steadily decreased. The effective
reproductive number however decreases until mid-April and starts to increase
again in most provinces over the summer. A seasonal increase in temperature is
present as the data was collected from late winter to mid-summer, while
humidity and PM2.5 do not follow any obvious trends. Descriptive statistics of
our data are available in Table 1.

| Mean | SD | Min | Max
---|---|---|---|---
Daily Cases | 42.29 | 76.93 | 0 | 868
Daily Tests | 968.93 | 1136.27 | 0 | 18256
Temperature (oC) | 18.06 | 6.60 | 0.50 | 36
Relative Humidity (%) | 63.84 | 17.67 | 11.50 | 100
PM2.5 ($\mu$g/m3) | 48.52 | 21.67 | 5 | 172
Mobility Decrease (%) | 15.57 | 12.31 | -9 | 47
Rt | 0.94 | 0.24 | 0.57 | 2.12
Table 1: Mean, standard deviation, minimum, and maximum of each variable
considered across all cities and days, including the estimated value of Rt.
Figure 2: The time series of daily cases, temperature, humidity, and PM2.5 in
23 Italian provinces. Each line represents one province and the colors are
grouped by region.

###  2.3 Statistical Model

Following other studies on the effects of ambient environmental variables on
viral transmission (air_pollution_and_mortality; dengue_gam;
temperature_covid_gam_brazil; temperature_covid_gam_china; leishmaniasis_gam;
feng2016impact), we used a generalized additive model (GAM) to analyze the
COVID-19 pandemic. The GAM is a semi-parametric extension of the generalized
linear model in which the response variable (relative transmission rate)
relies on a linear combination of smooth functions of the predictor variables
(humidity, PM2.5, temperature, and mobility). Our GAM was defined using a
Gaussian distribution and the following equation:

$\log(R_{t,i})=s(\text{Mobility}_{t,i})+s(\text{Temp}_{t,i})+s(\text{Hum}_{t,i})+s(\text{PM}_{2.5_{t,i}})+\lambda_{i},$
(1)

where $t$ denotes the date, $i$ denotes the province, and $s(\cdot)$ denotes
the thin plate basis spline used to smooth the data. Parameters of the splines
were estimated using the GCV method, and the splines were constrained to have
at most five degrees of freedom. This constraint was imposed to minimize
overfitting while accurately describing the trends of the data. The results of
the model were robust to changes in the degrees of freedom, signifying that we
effectively captured the trends in the data.

We chose our predictor variables to be the percent decrease in mobility
(Mobilityt,i), median temperature (Tempt,i), median humidity (Humt,i), and
median PM2.5 (PM${}_{2.5_{t,i}}$). We also included the fixed effect of each
province with the intercept $\lambda_{i}$ to account for local differences.
Finally, we selected our response variable to be the log-transformed effective
reproductive number ($R_{t,i}$). Choosing our response variable in this way
allowed us to avoid lagging any of our predictor variables, as the effective
reproductive number is already adjusted for the delay between the day of
infection and the day the case was reported (rtlive2020). Furthermore,
choosing the logarithm transformation allows the model to be interpreted more
naturally as we can see the effects of the predictor variables on a relative
scale rather than an absolute scale (GelmanHill).

The model was tested for robustness in two ways. First, we constructed the
model for each province individually in order to ensure that the trends
present in the global data set were also present in local data sets. Next, we
cross-validated the model on our data by removing one province from the global
data set, training the model on the 22 remaining provinces and using this
model to predict $R_{t}$ for the removed province. We then repeated this
process for each province in order to ensure that no province had an
overstated effect on the model.

The GAM was computed using the mgcv package in R (mgcv).

##  3 Results

Output from the GAM in Equation (1) is shown in Figure 3. These plots depict
the partial effect of each predictor variable on the effective reproductive
number $R_{t}$ determined by the evaluate$\\_$smooth function in the gratia
package in R (gratia). The 95% confidence intervals are shown in gray.
Negative values indicate a decrease in $R_{t}$, while positive values indicate
an increase in $R_{t}$. If the confidence interval includes zero, there is no
effect.

The percent change in mobility had a positive effect on $R_{t}$. This is a
very intuitive result; as people stay at home more, the effective reproductive
number decreases. Increased temperature showed a negative effect on the
effective reproductive number up until around 25$\degree$C. Humidity had a
somewhat positive effect when under 50$\%$ and a negligible effect at greater
values. PM2.5 had no effect at concentrations between 0 and 70 $\mu$g/m3, but
showed a strong positive effect for larger values. All predictor variables
were significant with p-values of less than $2\times 10^{-16}$. The model had
an adjusted coefficient of determination R${}^{2}=0.67$.

Figure 3: Each plot shows the partial effects of change in mobility, PM2.5,
temperature, and humidity on the expected value of $\log(R_{t})$. The 95%
confidence intervals are depicted in gray. $R_{t}$ increases as mobility
increases, PM2.5 has negligible effect for values between $0-70\mu$g/m3 but
displays a strong positive effect at higher values, temperature has a negative
relationship with $R_{t}$, and humidity has a slight positive effect under 50%
but a negligible effect at higher values.

The model was fit to each province individually which revealed similar
qualitative effects for each predictor as observed in the global data set
(Figure 4). Variability at higher PM2.5 values exists due to a small number of
observations per province. Firenze (Toscana), Napoli (Campagnia) and Roma
(Lazio) do not exhibit the same increasing trend for PM2.5 as the other
provinces but all had very few observations above 75 $\mu\text{g/m}^{3}$ and
no observations above 100 $\mu\text{g/m}^{3}$. Similarly, variability in the
effect of humidity below 50% can also be explained by a small sample size in
this range per province.

Figure 4: Each line shows the partial effect for each predictor variable when
fitting a GAM to each province separately. Lines were colored by region which
shows that the biggest outliers come from regions with only a small number of
observations.

Additionally, the model was trained on the data with one province removed at a
time and the mean squared error was calculated after predicting on the removed
province (Figure 5). Prediction error was lower in regions where we had data
from multiple provinces, and higher in regions where less data was available.
Illustrating this, Napoli, Roma, and Verona were the only provinces in their
respective regions from which we collected data, and they had the highest mean
squared prediction errors.

Figure 5: A geospatial heat map of mean squared prediction error from the
cross validation process. The error is lowest in the Lombardi region, where we
had the most observations, and highest in the regions of Campagni, Lazio, and
Veneto where we used data from one province each.

##  4 Discussion

The COVID-19 pandemic has affected millions of people around the world. Poor
air quality is also considered a public health issue, associated with long-
term health effects such as weakened cardiovascular health and decreasing life
expectancy (delfino2005potential; krewski2009evaluating;
lelieveld2015contribution; janssen2013short). Past studies have also shown an
association between high concentrations of particulate matter less than 2.5
$\mu$m in diameter and an increase in influenza transmission (su2019short).
Thus we hypothesize that even short-term exposure to elevated PM2.5 may lead
to a higher risk of COVID-19 transmission. Although any analysis to understand
transmission is confounded by many variables, such as temperature, humidity,
social distancing or the quality of data, our results show that air pollution
may be associated with a higher effective reproduction rate of the virus. We
showed the generalized additive model’s results to be robust. When each
province was removed from the data pool, a GAM constructed from the remaining
provinces could predict case counts in the excluded province.

Our study is novel in its use of effective reproductive number as the
dependent variable for the GAM. Although this parameter $R_{t}$ is commonly
calculated to estimate the growth of a pandemic (why_rt; why_rt2; why_rt3),
our use of it in a GAM analyzing the impact of other variables on disease
spread is novel. As we were interested in a daily change in virus
transmission, we chose to use the assumptions of an established model to
account for the many unknowns in the data, such as the lag between contracting
the virus and showing symptoms, and the one between taking a test and
receiving a positive result. Furthermore, we use a longer timescale than
comparable models of the COVID-19 pandemic, including both the initial wave of
cases as well as the subsequent decline. Finally, we account for changes in
human behavior during the pandemic using Google mobility data, while the
majority of contemporary papers have not incorporated this or any similar
metric.

We found statistically significant relationships of mobility, temperature,
humidity, and PM2.5 on COVID-19 transmission. It is well established that the
reduction in mobility as a result of national lockdowns has been very
effective at reducing the spread of COVID-19 (effect_of_mobility;
china_mobility_aqi_gam), which is clear from our results as well. The
magnitude of this effect explains why it was essential to include in our
analysis. We also found an increase in temperature was associated with a
decrease in $R_{t}$. We are hesitant to ascribe excessive significance to this
finding given the strong serial correlation of temperature and the limited
data we have available. Nevertheless, we can still compare our findings to
those of other studies of COVID-19. 166_countries used a linear model on data
from 166 countries to note a negative correlation between temperature and
COVID-19 cases. They also found a significant correlation between humidity and
COVID-19 cases, albeit with a noticeably smaller effect than that of
temperature. Both of these are supported by our results. Furthermore,
temperature_covid_gam_brazil used a non-linear model in Brazil to find a
noticeable decrease in transmission as temperature increased from
$17^{\circ}$C to $24^{\circ}$C and little change between $24-27^{\circ}$C,
which was similar to our findings. Studies of China such as those by
temperature_covid_gam_china and temperature_covid_gam_china2 have found the
opposite result. They report positive nonlinear trends between temperature and
disease incidence in January and February. Possible reasons for this
disagreement may be model selection, confounding variables, or a limited range
of data (166_countries).

Our results regarding PM2.5 were also largely consistent with previous studies
on COVID-19 and other respiratory illnesses. feng2016impact, liang2014pm,
janssen2013short, su2019short, jiang2020effect, yongjian2020association, and
china_mobility_aqi_gam observed that PM2.5 was positively correlated to the
transmission of various viruses including coronaviruses and influenza viruses.
Although most of these studies used linear models to establish a correlation,
feng2016impact reported a nonlinear trend with little correlation between
infection rates for PM2.5 values below 70 $\mu$g/m3, followed by a linearly
increasing trend for PM2.5 values between 70 and 150 $\mu$g/m3, a similar
phenomenon to the one we observed. Overall, it seems the majority of studies
agree that a statistically significant relationship between PM2.5 and viral
transmission exists. The common hypothesis for this is two-fold. First, fine
particulate matter has been reported to damage the respiratory system,
possibly resulting in “a temporary immunosuppressive pulmonary
microenvironment" (feng2016impact; inhalation_toxicology;
inhalation_toxicology2; inhalation_toxicology3; inhalation_toxicology4).
Second, PM2.5 is sufficiently small that it can stay suspended in the air for
prolonged stretches of time, making it, and any pathogen that may have bound
itself to it, susceptible to inhalation. In fact, a recent finding from the
Italian province of Bergamo provides evidence that SARS-CoV-2 does bind to
particulate matter in the air (bergamo_evidence). Both of these hypotheses
could imply that particulate matter allows greater exposure to the virus as
well as a greater susceptibility to infection following exposure. However,
more research is needed to understand the underlying mechanisms and the degree
of additional risk.

As pollution levels rise around the world (shaddick2020half), it is important
to study their impact on respiratory disease. Even now, events such as
wildfires in the western United States drive PM2.5 levels to extreme highs
(liu2016particulate; xie2020summer). Thus, our findings may predict an
association between the fires and higher incidences of COVID-19 transmission.
Furthermore, as the climate continues to change, this association may be
important to note in future outbreaks of respiratory illnesses.

Another pervasive issue in the pandemic is the disproportional incidence of
COVID-19 on low-income and ethnic communities, particularly in the United
States (racial_disparities; racial_disparities2; race_covid_air).
Understanding the causes of this disparity would allow for policymakers to
implement more effective measures of controlling the spread in these
communities. To this end, it has been established that low-income and ethnic
communities are disproportionately exposed to high levels of PM2.5
(race_covid_air; race_air). The results of our study indicate that it is
possible that some of the increased disease transmission in these communities
may be attributed to higher levels of air pollution. Unfortunately due to a
lack of available data, we were unable to directly investigate this using a
subgroup analysis; however, we bring it up as an important issue that warrants
further consideration.

In addition to PM2.5, temperature, and humidity, it is possible to consider
other environmental variables. Other studies have identified possible
relationships between COVID-19 transmission and PM10, SO2, NO2, and O3
(china_mobility_aqi_gam; jiang2020effect; Fronza2020). We did not include
these in our study for a variety of reasons. There were very few days where
any province reported unhealthy NO2 levels ($>$150 ppb), and as a result, we
unsurprisingly found no strong effect. Likewise, O3 measurements were missing
from a significant number of provinces, and the O3 that was reported seldom
reached unhealthy ($>$100 ppb) levels. PM10 and SO2 measurements were
similarly unavailable for a large number of provinces and thus were omitted
from study.

There are several limitations to our study. While we attempt to control for
the seasonal variation of atmospheric variables in our model, an entire
seasonal cycle of the COVID-19 pandemic has not been observed, making a full
seasonal analysis infeasible at this time. Many confounding variables may
influence PM2.5 levels, and they are difficult to disentangle. Furthermore,
since our study focuses on Italy, we do not know if these findings would be
consistent in other areas. Future work would include expanding this study to
multiple countries in different regions, such as the southern hemisphere, and
to countries with a less homogeneous government response to test if our
results hold on a global scale.

##  5 Conclusion

Our results suggest that mobility, temperature, humidity, and short-term
exposure to PM2.5 are all possible risk factors for COVID-19. While the
effects of temperature, humidity, and mobility have all been subjects of
large-scale global studies on the pandemic (effect_of_mobility;
166_countries), PM2.5 has received less attention. However, because PM2.5 can
change significantly due to both long-term, predictable events (such as an
increase in PM2.5 in winter months) as well as short-term, less predictable
events (such as wildfires), it is important to try to understand how it may
impact the pandemic. We find that PM2.5 has a negligible correlation with
$R_{t}$ for healthy to moderate air quality, but as PM levels increase to
unhealthy amounts, a strong positive correlation emerges. This provides
evidence for the long-standing hypothesis that short-term exposure to air
pollution is a risk factor for respiratory illnesses. Although the results of
our study are consistent with other works, this should still be considered a
preliminary study as larger scale data sets are needed to confirm a global
trend.

## Supplementary Materials

The data and R code used to generate these results, as well as figures for
individual cities in the validation process, are included in the supplementary
files. A README is provided to describe the data available, how to generate
each figure presented in the paper, and where important variables should be
found in the code.

## Acknowledgements

We thank Dr. Richard McGehee, Dr. Mary Silber, and Dr. Mary Lou Zeeman for
their contributions advising this effort as well as Dr. Jack O’Brien and Dr.
Douglas Nychka for discussions regarding the development of the statistical
model. This work was supported by the NSF, the American Institute of
Mathematics, and the Mathematics and Climate Research Network.

## References

  * Brandt et al. (2020) Brandt EB, Beck AF, Mersha TB (2020).  Air pollution, racial disparitiesm and COVID-19 mortality.  _J Allergy Clin Immunol_ , 146(1): 61–63.
  * Carslaw (2020) Carslaw D (2020).  _Worldmet: Import Surface Meteorological Data from NOAA Integrated Surface Database (ISD)_.  R package version 0.9.0.
  * Ciencewicki and Jaspers (2007) Ciencewicki J, Jaspers I (2007).  Air pollution and respiratory viral infection.  _Inhalation Toxicology_ , 19(14): 1135–1146.
  * Colon-Gonzalez et al. (2013) Colon-Gonzalez FJ, Fezzi C, Lake IR, Hunter PR (2013).  The effects of weather and climate change on dengue.  _PLoS Neglected Tropical Diseases_ , 7: e2503.
  * Contini and Costabile (2020) Contini D, Costabile F (2020).  Does air pollution influence COVID-19 outbreaks?
  * Delfino et al. (2005) Delfino RJ, Sioutas C, Malik S (2005).  Potential role of ultrafine particles in associations between airborne particle mass and cardiovascular health.  _Environmental Health Perspectives_ , 113(8): 934–946.
  * Donaldson and Tran (2002) Donaldson K, Tran CL (2002).  Inflammation caused by particles and fibers.  _Inhalation Toxicology_ , 14: 5–27.
  * Feng et al. (2016) Feng C, Li J, Sun W, Zhang Y, Wang Q (2016).  Impact of ambient fine particulate matter (PM 2.5) exposure on the risk of influenza-like-illness: a time-series analysis in Beijing, China.  _Environmental Health_ , 15(1): 17.
  * Fronza et al. (2020) Fronza R, Lusic M, Schmidt M, Lucic B (2020).  Spatial–temporal variations in atmospheric factors contribute to SARS-CoV-2 outbreak.  _Viruses_ , 12(6): 588.
  * Gelman and Hill (2007) Gelman A, Hill J (2007).  _Data analysis using regression and multilevel/hierarchical models_ , volume Analytical methods for social research.  Cambridge University Press, New York.
  * Hooper et al. (2020) Hooper MW, Nápoles AM, Pérez-Stable EJ (2020).  COVID-19 and racial/ethnic disparities.  _JAMA_ , 323(24): 2466–2467.
  * Jaffe et al. (2008) Jaffe D, Hafner W, Chand D, Westerling A, Spracklen D (2008).  Interannual variations in PM2.5 due to wildfires in the western United States.  _Environmental Science & Technology_, 42(8): 2812–2818.
  * Janssen et al. (2013) Janssen N, Fischer P, Marra M, Ameling C, Cassee F (2013).  Short-term effects of PM2. 5, PM10 and PM2.5–10 on daily mortality in the Netherlands.  _Science of the Total Environment_ , 463: 20–26.
  * Jiang et al. (2020) Jiang Y, Wu XJ, Guan YJ (2020).  Effect of ambient air pollutants and meteorological variables on COVID-19 incidence.  _Infection Control & Hospital Epidemiology_, 41(9): 1011–1015.
  * Kim and Kang (1995) Kim CS, Kang TC (1995).  Comparative measurement of lung deposition of inhaled fine particles in normal subjects and patients with obstructive airway disease.  _American Journal of Respiratory and Critical Care Medicine_ , 155: 899––905.
  * Kraemer et al. (2020) Kraemer MU, Yang CH, Gutierrez B, Wu CH, Klein B, Pigott DM, et al. (2020).  The effect of human mobility and control measures on the COVID-19 epidemic in China.  _Science_ , 368: 493–497.
  * Krewski (2009) Krewski D (2009).  Evaluating the effects of ambient air pollution on life expectancy.
  * Laden et al. (2000) Laden F, Neas LM, Dockery DW, Schwartz J (2000).  Association of fine particulate matter from different sources with daily mortality in six US cities.  _Environmental Health Perspectives_ , 108(10): 941–947.
  * Lauer et al. (2020) Lauer SA, Grantz KH, Bi Q, Jones FK, Zheng Q, Meridith HR, et al. (2020).  The incubation period of coronavirus disease 2019 (COVID-19) from publicly reported confirmed cases: Estimation and application.  _Annals of Internal Medicine_ , 172(9): 577–582.
  * Lee et al. (2014) Lee GI, Saravia J, You D, Shrestha B, Jaligama S, Hebert VY, et al. (2014).  Exposure to combustion generated environmentally persistent free radicals enhances severity of influenza virus infection.  _Particle and Fibre Toxicology_ , 11: 11––57.
  * Lee et al. (2020) Lee GI, Saravia J, You D, Shrestha B, Jaligama S, Hebert VY, et al. (2020).  Effects of temperature and humidity on the daily new cases and new deaths of COVID-19 in 166 countries.  _Science of the Total Environment_ , 729: 139051.
  * Lelieveld et al. (2015) Lelieveld J, Evans JS, Fnais M, Giannadaki D, Pozzer A (2015).  The contribution of outdoor air pollution sources to premature mortality on a global scale.  _Nature_ , 525(7569): 367–371.
  * Li et al. (1996) Li XY, Gilmour PS, Donaldson K, MacNee W (1996).  Free radical activity and pro-inflammatory effects of particulate air pollution (PM10) in vivo and in vitro.  _Thorax_ , 51: 1216–1222.
  * Liang et al. (2014) Liang Y, Fang L, Pan H, Zhang K, Kan H, Brook JR, et al. (2014).  PM 2.5 incBeijing–temporal pattern and its association with influenza.  _Environmental Health_ , 13(1): 102.
  * Lipsitch et al. (2011) Lipsitch M, Finelli L, Herrenan RT, Leung GM, Redd SC (2011).  Improving the evidence base for decision making during a pandemic: The example of 2009 Influenza A/H1N1.  _Biosecurity and Bioterrorism: Biodefense Strategy, Practice, and Science_ , 9: 89–115.
  * Liu et al. (2016) Liu JC, Mickley LJ, Sulprizio MP, Dominici F, Yue X, Ebisu K, et al. (2016).  Particulate air pollution from wildfires in the western US under climate change.  _Climatic change_ , 138(3-4): 655–666.
  * Ma et al. (2020) Ma Y, Zhao Y, Liu J, He X, Wang B, Fu S, et al. (2020).  Effects of temperature variation and humidity on the death of COVID-19 in Wuhan, China.  _Science of the Total Environment_ , 724: 138226.
  * Nishiura et al. (2020) Nishiura H, Linton NM, Akhmetzhanov AR (2020).  Serial interval of novel coronavirus (COVID-19) infections.  _International Journal of Infectious Diseases_ , 93: 284–286.
  * O’Dell et al. (2019) O’Dell K, Ford B, Fischer EV, Pierce JR (2019).  Contribution of wildland-fire smoke to US PM2. 5 and its influence on recent trends.  _Environmental Science & technology_, 53(4): 1797–1804.
  * Pan et al. (2020) Pan A, Liu L, Wang C, Guo H, Hao X, Wang Q, et al. (2020).  Association of public health interventions with the epidemiology of the COVID-19 outbreak in Wuhan, China.  _JAMA Network_ , 323: 1915–1923.
  * Peng et al. (2005) Peng RD, Dominici F, Louis TA (2005).  Model choice in time series studies of air pollution and mortality.  _Journal of the Royal Statistical Society_ , 169: 179–203.
  * Prata et al. (2020) Prata DN, Rodrigues W, Bermejo PH (2020).  Temperature significantly changes COVID-19 transmission in (sub)tropical cities of Brazil.  _The Science of the Total Environment_ , 729: 138862.
  * Richardson et al. (2020) Richardson S, Hirsch JS, Narasimhan M, Crawford JM, McGinn T, Davidson KW, et al. (2020).  Presenting characteristics, comorbidities, and outcomes among 5700 patients hospitalized with COVID-19 in the New York City area.  _JAMA_ , 323(20): 2052–2059.
  * Rubin et al. (2020) Rubin D, Huang J, Fisher BT, Gasparrini A, Tam V, Song L, et al. (2020).  Association of social distancing, population density, and temperature with the instantaneous reproduction number of SARS-CoV-2 in counties across the United States.  _JAMA Network Open_ , 3: e2016099.
  * Schwartz et al. (1996) Schwartz J, Dockery DW, Neas LM (1996).  Is daily mortality associated specifically with fine particles?  _Journal of the Air & Waste Management Association_, 46(10): 927–939.
  * Setti et al. (2020) Setti L, Passarini F, De Gennaro G, Barbieri P, Perrone MG, Borelli M, et al. (2020).  SARS-CoV-2RNA found on particulate matter of Bergamo in northern Italy: First evidence.  _Environmental Research_ , 188: 109754.
  * Shaddick et al. (2020) Shaddick G, Thomas M, Mudu P, Ruggeri G, Gumy S (2020).  Half the world’s population are exposed to increasing air pollution.  _NPJ Climate and Atmospheric Science_ , 3(1): 1–5.
  * Simpson (2021) Simpson G (2021).  _Graceful ’ggplot’-Based Graphics and Other Functions for GAMs Fitted Using ’mgcv’_.  R package version 0.5.1.
  * Su et al. (2019) Su W, Wu X, Geng X, Zhao X, Liu Q, Liu T (2019).  The short-term effects of air pollutants on influenza-like illness in Jinan, China.  _BioMed Central Public Health_ , 19(1): 1–12.
  * Systrom et al. (2020) Systrom K, Vladek T, Krieger M (2020).  Rt.live.  _GitHub repository_.
  * Talmoudi et al. (2017) Talmoudi K, Bellali H, Ben-Alaya N, Saez M, Malouche D, Chahed MK (2017).  Modeling zoonotic cutaneous leishmaniasis incidence in central Tunisia from 2009-2015: Forecasting models using climate variables as predictors.  _PLoS Neglected Tropical Diseases_ , 11: e0005844.
  * Tessum et al. (2019) Tessum CW, Apte JS, Goodkind AL, Muller NZ, Mullins KA, Paolella DA, et al. (2019).  Inequality in consumption of goods and services adds to racial-ethnic disparities in air pollution exposure.  _PNAS_ , 116(13): 6001–6006.
  * Tucker (2000) Tucker WG (2000).  An overview of PM2.5 sources and control strategies.  _Fuel Processing Technology_ , 65: 379–392.
  * Wood (2011) Wood SN (2011).  Fast stable restricted maximum likelihood and marginal likelihood estimation of semiparametric generalized linear models.  _Journal of the Royal Statistical Society (B)_ , 73(1): 3–36.
  * Xie and Zhu (2020) Xie J, Zhu Y (2020).  Association between ambient temperature and COVID-19 infection in 122 cities from China.  _The Science of the Total Environment_ , 724: 138201.
  * Xie et al. (2020) Xie Y, Lin M, Horowitz LW (2020).  Summer PM2. 5 pollution extremes caused by wildfires over the western United States during 2017–2018.  _Geophysical Research Letters_ , 47(16): e2020GL089429.
  * Xing et al. (2016) Xing YF, Xu YH, Shi MH, Lian YX (2016).  The impact of PM2. 5 on the human respiratory system.  _Journal of Thoracic Disease_ , 8(1): E69.
  * Xu et al. (2020) Xu B, Gutierrez B, Mekaru S, Sewalk K, Goodwin L, Loskill A, et al. (2020).  Epidemiological data from the COVID-19 outbreak, real-time case information.  _Scientific Data_ , 7.
  * Yancy (2020) Yancy CW (2020).  COVID-19 and African Americans.  _JAMA_ , 323(19): 1891–1892.
  * Yang et al. (2020) Yang J, Zheng Y, Gou X, Pu K, Chen Z, Guo Q, et al. (2020).  Prevalence of comorbidities and its effects in patients infected with SARS-CoV-2: a systematic review and meta-analysis.  _International Journal of Infectious Diseases_ , 94: 91–95.
  * Zhu et al. (2020a) Zhu N, Zhang D, Wang W, Li X, Yang B, Song J, et al. (2020a).  A novel coronavirus from patients with pneumonia in china, 2019.  _New England Journal of Medicine_ , 382: 727–733.
  * Zhu et al. (2020b) Zhu Y, Xie J, Huang F, Cao L (2020b).  Association between short-term exposure to air pollution and covid-19 infection: Evidence from China.  _Science of the Total Environment_ , 138704.
  * Zhu et al. (2020c) Zhu Y, Xie J, Huang F, Cao L (2020c).  The mediating effect of air quality on the association between human mobility and COVID-19 infection in China.  _Environmental Research_ , 189: 109911.


# Sparsity-Based Kalman Filters for Data Assimilation††thanks: This work was
supported in part by U.S. Naval Research Laboratory - Monterey, CA

Wei Kang and Liang Xu Wei Kang is with Faculty of Applied Mathematics, Naval
Postgraduate School, Monterey, CA, USA<EMAIL_ADDRESS>Xu is with Naval
Research Laboratory, Monterey, CA, USA<EMAIL_ADDRESS>
###### Abstract

Several variations of the Kalman filter algorithm, such as the extended Kalman
filter (EKF) and the unscented Kalman filter (UKF), are widely used in science
and engineering applications. In this paper, we introduce two algorithms of
sparsity-based Kalman filters, namely the sparse UKF and the progressive EKF.
The filters are designed specifically for problems with very high dimensions.
Different from various types of ensemble Kalman filters (EnKFs) in which the
error covariance is approximated using a set of dense ensemble vectors, the
algorithms developed in this paper are based on sparse matrix approximations
of error covariance. The new algorithms enjoy several advantages. The error
covariance has full rank without being limited by a set of ensembles. In
addition to the estimated states, the algorithms provide updated error
covariance for the next assimilation cycle. The sparsity of error covariance
significantly reduces the required memory size for the numerical computation.
In addition, the granularity of the sparse error covariance can be adjusted to
optimize the parallelization of the algorithms.

##  1 INTRODUCTION

For dynamical systems with numerical models, data assimilation is an
estimation process of combining observational data with a numerical model to
obtain an estimate of the system’s state. Data assimilation is essential to
numerical weather prediction (NWP). The state estimate can be considered as an
interpolation of the sparse observational data; and it is used as the initial
condition for the numerical forecast process. If the dimension is relatively
low and the data set is small, various linear and nonlinear estimators can be
found in the literature that have optimal or suboptimal performances. However,
to assimilate big data sets with models that have high dimensions, such as
those in operational NWP systems with tens of millions of variables, achieving
reliable state estimates and error probability distributions is a challenging
problem that have been studied for decades with a huge literature.

There are two categories of methods widely used in NWP, namely variational
methods and the ensemble Kalman filter (EnKF) [8, 3]. The former is based on a
weighted least-square optimization, such as the four dimensional variational
data assimilation (4D-Var) in a fixed time window or the three dimensional
version (3D-Var) that excludes the time variable. The EnKF algorithm is based
on the Kalman filter except that the error covariance is approximated using a
set of state ensembles. 4D-Var methods are used in operational NWP systems by
many meteorological centers. While it serves as an effective method of data
assimilation, 4D-Var algorithms have difficulty to explicitly track the
evolution of error covariance within its estimation process due to outrageous
computational costs and input/output (I/O) loads required by matrices of
extremely high dimensions. EnKF, on the other hand, updates information about
the error covariance in the form of ensembles. However, the number of ensemble
states is significantly smaller than the number of state variables. As a
result, the rank deficiency of error covariance tends to deteriorate the
integrity of the estimation process unless remedies to the algorithm, such as
localization and covariance inflation, are applied.

Different types of Kalman filters have been developed and widely used in
science and engineering applications, such as the EnKF, the extended Kalman
filter (EKF) and the unscented Kalman filter (UKF). In this paper, we
introduce two algorithms of sparsity-based Kalman fillters, namely the sparse
UKF and the progressive EKF. The goal of the work is to explore innovative
ideas that take the advantage of the (almost) sparsity structure of matrices
so that analysis and error covariance can be updated effectively and
efficiently without the drawback of rank deficiency. The granularity of
subproblems for the purpose of algorithm parallelization is also emphasized in
the method. The filters are developed specifically for problems with very high
dimensions. Different from EnKFs in which the error covariance is represented
using a set of dense ensemble vectors, the new algorithms in this paper are
based on a sparse but full rank matrix as an approximation of the error
covariance. This is made possible because of two assumptions: (a) the error
covariance is approximately a sparse matrix; (b) the system model is component
based, i.e. the state vectors are divided into components that can be computed
independently in parallel.

##  2 Sparse UKF

Consider a dynamical system model in which the state variable is $x(t)$, where
$t=1,2,3,\cdots$ represents time steps. The value of observation at $t=k$ is
denoted by $y(k)$. The system model is defined as follows,

$\begin{array}[]{lllllllll}\begin{aligned} x(k)&={\cal
M}(x(k-1))+\eta_{k-1},&&x(k),\eta_{k-1}\in\mathbb{R}^{n},\\\ y(k)&={\cal
H}(x(k))+\delta_{k},&&y_{k},\delta_{k}\in\mathbb{R}^{m},\end{aligned}\end{array}$
(1)

where $\eta_{k-1}$ is a random variable representing the model error. Its
covariance is $Q$. The observational error, $\delta_{k}$, has a covariance
$R$. In data assimilation, the goal is to estimate the value of $x(k)$ given
the observations $y(1),y(2),\cdots,y(k)$ and the model (1). If (1) is linear
and if all random variables are Gaussian, then the Kalman filter is an optimal
state estimator. For nonlinear systems with non-Gaussian random errors,
various types of Kalman filters exist in the literature with successful
applications in science and engineering. If a system has a very high
dimension, the conventional form of Kalman filter based on a dense error
covariance is not applicable. In this section, we introduce an algorithm that
is a variation of UKF for problems with approximately sparse error
covariances.

In a sparse matrix/vector, most entries are zeros. In this paper, we use an
underbar to emphasize that a vector or matrix is sparse, for instance
$\underline{P}$ and $\underline{x}$. A sparse vector is associated with an
index set, denoted by ${\cal I}$, consisting of the indices of nonzero
entries. For a sparse matrix, columns may have different numbers of nonzero
entries. The largest such number is denoted by $N_{sp}$. In sparsity-based
algorithms, a full model evaluation is not always necessary. Using a
component-based model can significantly reduce the computational load. In the
notation, a component-based model has three inputs: the sparse state variable,
its index, and the index of the output state. More specifically,

$\begin{array}[]{lllllllll}\underline{x}(k)={\cal M}(\underline{x}(k-1);{\cal
I}_{1};{\cal I}_{2}),\\\ \end{array}$ (2)

where ${\cal I}_{1}$ is the index set of the sparse vector
$\underline{x}(k-1)$ and ${\cal I}_{2}$ is the index set of
$\underline{x}(k)$. The model evaluates only the entries with indices in
${\cal I}_{2}$, setting all other entries as zeros. Note that
$\underline{x}(k)$ is different from the full state variable $x(k)$ because
the later is, in general, a dense vector with mostly nonzero entries.
Therefore, it is important to specify the index set ${\cal I}_{2}$ of the
sparse vector $\underline{x}(k)$ to be evaluated using a component-based
model. For simplicity, we often omit ${\cal I}_{1}$ in the notation, i.e.

$\begin{array}[]{lllllllll}\underline{x}(k)={\cal M}(\underline{x}(k-1);{\cal
I}),\end{array}$ (3)

where ${\cal I}$ is the same as ${\cal I}_{2}$ in (2). Algebraic operations
between sparse vectors, such as addition and dot product, are defined in the
same way as dense vectors. Thus, we may conduct operations between sparse
vectors and dense vectors, such as adding a sparse vector to a dense vector
$\underline{x}+y$ as long as both vectors have the same dimension. A new
operation, called merging, between a sparse vector and a dense vector is
defined as follows,

$\begin{array}[]{lllllllll}z=\underline{x}\triangleright
y,&\left\\{\begin{array}[]{lll}i\mbox{th component of }z=i\mbox{th component
of }x,&\mbox{if }i\in{\cal I}.\\\ i\mbox{th component of }z=i\mbox{th
component of }y,&\mbox{if }i\not\in{\cal I}.\end{array}\right.\end{array}$ (4)

If an operation has an underbar, it means that the evaluation is carried out
only at a given index set. For instance, given two sparse matrices
$\underline{A}$ and $\underline{B}$, then $\underline{A}*\underline{B}$ is a
different matrix from $\underline{\underline{A}*\underline{B}}$. The former is
the conventional matrix multiplication between two sparse matrices; the later
is a matrix multiplication in which the entries in a given index set are
evaluated and all other entries are set to be zeros. Other operations, such as
$\underline{\sqrt{\underline{P}}}$, are defined similarly. A summary of
notations is summarized in the following table.

Notation | Definition | Notation | Definition
---|---|---|---
$x$ | state variable | $y$ | observation variable
${\cal M}$ | model function | ${\cal H}$ | observation operator
$n$ | state space dimension | $t=1,2...$ | (discrete) time variable
$x^{\sigma}_{i}$ |  $\sigma$-point at $t=k-1$ |  |
$x^{b}_{i}$ | background - state vector | $y^{b}_{i}$ | output of observation operator ${\cal H}(x^{b}_{i})$
$\bar{x}^{b}$ | average of $x^{b}_{i}$ | $\bar{y}^{b}$ | average of $y^{b}_{i}$
$P^{b}$ | background - error covariance |  |
$x^{a}$ | analysis - state vector | $P^{a}$ | analysis - error covariance
Table 1: Notations

###  2.1 The UKF

The unscented Kalman filter has been increasingly popular in engineering
applications since its introduction about twenty years ago [4, 5]. In a UKF
algorithm, the error covariance is propagated with the dynamics using a set of
vectors, or $\sigma$-points denoted by $x^{\sigma}$. Their definition is given
in (5)-(6). The $\sigma$-points are computed at each time step using a square
root of the error covariance. In most UKF applications, $\sigma$-points are
computed using either Cholesky factorization or matrix diagonalizations. In
the notation, a variable with a superscript ’$a$’, such as $x^{a}$, represents
the analysis value of the variable, i.e., the updated value based on
observations. A variable with a superscript ’$b$’, such as $y^{b}$, represents
the background, i.e., the propagated value of analysis using the system model.
The algorithm is summarized as follows. At $t=k-1$, suppose we have the
analysis and error covariance as well as its square root

$\begin{array}[]{lllllllll}x^{a}(k-1),\;P^{a}(k-1),\\\
X^{a}(k-1)=\sqrt{(n+\kappa)P^{a}(k-1)},\end{array}$ (5)

where $\kappa$ is a scaling factor for the fine tuning of the higher order
moments of the approximation error [4]. How to tune the value of $\kappa$ for
a sparsity-based UKF is an open problem that needs further study. In this
paper, $\kappa=0$ is used in all examples. A set of $\sigma$-points is
generated as follows,

$\begin{array}[]{lllllllll}x^{\sigma}_{0}(k-1)=x^{a}(k-1),\\\
x^{\sigma}_{i}(k-1)=x^{a}(k-1)+X^{a}_{i}(k-1),&1\leq i\leq n,\\\
x^{\sigma}_{i}(k-1)=x^{a}(k-1)-X^{a}_{i}(k-1),&n+1\leq i\leq 2n.\\\
\end{array}$ (6)

The next step is to propagate the $\sigma$-points, which represent the
background at $t=k$. For simplicity, the index ’$k$’ of all variables in the
$k$th time-step is omitted.

$\displaystyle x^{b}_{i}$ $\displaystyle={\cal M}(x^{\sigma}_{i}(k-1)),$
$\displaystyle y^{b}_{i}$ $\displaystyle={\cal H}(x^{b}_{i}),$ $\displaystyle
0\leq i\leq 2n,$ (7) $\displaystyle\bar{x}^{b}$
$\displaystyle=\sum_{i=0}^{2n}w_{i}x^{b}_{i},$ $\displaystyle\bar{y}^{b}$
$\displaystyle=\sum_{i=0}^{2n}w_{i}y^{b}_{i},$

where the weights are defined as follows

$\begin{array}[]{lllllllll}w_{0}=\displaystyle\frac{\kappa}{n+\kappa},\;w_{i}=\frac{1}{2(n+\kappa)},\end{array}$
(8)

for $i=1,2,\cdots,2n$. Define the variations

$\begin{array}[]{lllllllll}X^{b}_{i}=x^{b}_{i}-\bar{x}^{b},\;Y^{b}_{i}=y^{b}_{i}-\bar{y}^{b}.\end{array}$
(9)

The background covariances are

$\begin{array}[]{lllllllll}P^{b}=\displaystyle\sum_{i=0}^{2n}w_{i}X^{b}_{i}(X^{b}_{i})^{T}+Q,\\\
P_{xy}=\displaystyle\sum_{i=0}^{2n}w_{i}X^{b}_{i}(Y^{b}_{i})^{T},\\\
P_{yy}=\displaystyle\sum_{i=0}^{2n}w_{i}Y^{b}_{i}(Y^{b}_{i})^{T}+R.\end{array}$
(10)

The Kalman gain, $K$, satisfies the following equation,

$\begin{array}[]{lllllllll}KP_{yy}=P_{xy}.\end{array}$ (11)

The analysis is updated as follows

$\begin{array}[]{lllllllll}\begin{aligned}
x^{a}&=\bar{x}^{b}+K(y_{o}-\bar{y}^{b}),\\\
P^{a}&=P^{b}-K(P_{xy})^{T}.\end{aligned}\end{array}$ (12)

where $y_{o}$ is the observation at $t=k$. This completes one iteration of the
filter. For the next step, $t=k+1$, go back to (5) replacing the analysis by
the updated value of $x^{a}$ and $P^{a}$.

###  2.2 Sparse UKF

The square root factorization of a matrix is not unique. For large and sparse
matrices, various algorithms and their implementations on different computing
platforms have been studied for many years. The literature can be traced back
to the early days of electronic computers [2]. In the case of Cholesky
factorization, the square root of a sparse matrix is still sparse, although
the computation may require larger memory than the original matrix [1, 7].

A dense error covariance is intractable in computation for global models used
in NWP. In the proposed approach, we assume that $P$ and $\sqrt{P}$ are
approximately sparse. In the algorithm, they are replaced by their sparse
approximations, $\underline{P}$ and $\sqrt{\underline{P}}$. The indices of
nonzero entries in the $i$th column are denoted by ${\cal I}_{i}$ and ${\cal
I}^{\sigma}_{i}$, respectively. The sparsity index set of the forecast of
$\sigma$-points, i.e. the background, are denoted by ${\cal I}^{b}$. Becasue
$\underline{P}$ and $\sqrt{\underline{P}}$ are approximations of $P$ and
$\sqrt{P}$, the sparsity patterns do not have to be exact. In fact, in the
example of Lorenz-96 model presented in the next section, we assume that
${\cal I}_{i}$, ${\cal I}^{\sigma}_{i}$ and ${\cal I}^{b}$ equal to each other
although $\sqrt{\underline{P}}$ may have a different sparsity pattern from
that of $\underline{P}$.

Algorithm I (Sparse UKF)
Given the initial analysis,

$\begin{array}[]{lllllllll}x^{a}(k-1),\;\underline{P}^{a}(k-1).\end{array}$
(13)

Step 1. $\sigma$-points and forecast

$\underline{X}^{a}(k-1)=\underline{\sqrt{(n+\kappa)\underline{P}^{a}(k-1)}},\mbox{
sparsity index set }{\cal I}^{\sigma}$ (14)

and

$\displaystyle x^{b}_{0}$ $\displaystyle={\cal M}(x^{a}(k-1)),$ $\displaystyle
y^{b}_{0}$ $\displaystyle={\cal H}(x^{b}_{0}),$ (15)
$\displaystyle\underline{x}^{b}_{i}$ $\displaystyle={\cal
M}(x^{a}(k-1)+\underline{X}^{a}_{i}(k-1);{\cal I}^{b}_{i}),$ $\displaystyle
y^{b}_{i}$ $\displaystyle={\cal H}(\underline{x}^{b}_{i}\triangleright
x^{b}_{0}),$ $\displaystyle 1\leq i\leq n.$
$\displaystyle\underline{x}^{b}_{i}$ $\displaystyle={\cal
M}(x^{a}(k-1)-\underline{X}^{a}_{i}(k-1);{\cal I}^{b}_{i}),$ $\displaystyle
y^{b}_{i}$ $\displaystyle={\cal H}(\underline{x}^{b}_{i}\triangleright
x^{b}_{0}),$ $\displaystyle n+1\leq i\leq 2n.$ $\displaystyle\bar{x}^{b}$
$\displaystyle=\sum_{i=0}^{2n}w_{i}(\underline{x}^{b}_{i}\triangleright
x^{b}_{0}),$ $\displaystyle\bar{y}^{b}$
$\displaystyle=\sum_{i=0}^{2n}w_{i}y^{b}_{i}$

Step 2. Background covariances

$\begin{array}[]{lllllllll}\begin{aligned}
\underline{P}^{b}&=\displaystyle\sum_{i=0}^{2n}w_{i}\underline{(\underline{x}^{b}_{i}\triangleright
x^{b}_{0}-\bar{x})(\underline{x}^{b}_{i}\triangleright
x^{b}_{0}-\bar{x})^{T}}+\underline{Q},&\mbox{ sparsity index set }{\cal I},\\\
P_{xy}&=\displaystyle\sum_{i=0}^{2n}w_{i}(\underline{x}^{b}_{i}\triangleright
x^{b}_{0}-\bar{x})(y^{b}_{i}-\bar{y})^{T},\\\
P_{yy}&=\displaystyle\sum_{i=0}^{2n}w_{i}(y^{b}_{i}-\bar{y})(y^{b}_{i}-\bar{y})^{T}+R.\end{aligned}\end{array}$
(16)

Step 3. Kalman gain and analysis

$\begin{array}[]{lllllllll}\begin{aligned} KP_{yy}&=P_{xy},\\\
x^{a}&=\bar{x}^{b}+K(y_{o}-\bar{y}^{b}),\\\
\underline{P}^{a}&=\underline{P}^{b}-\underline{K(P_{xy})}^{T}+\gamma
I,&\mbox{ sparsity index set }{\cal I}.\end{aligned}\end{array}$ (17)

In (16), we assume $\underline{x}^{b}_{0}=x^{b}_{0}$. The constant term
$\gamma I$ in (17) is a diagonal matrix. The value of $\gamma$ is selected so
that $P^{a}$ is positive definite. The positive definiteness is guaranteed if
$\gamma$ is larger than $|\lambda_{\min}|$, where $|\lambda_{\min}|$ is the
smallest negative eigenvalue of

$\begin{array}[]{lllllllll}\underline{P}^{b}-\underline{K(P_{xy})}^{T}.\end{array}$
(18)

If the updated covariance matrix is positive definite, then $\gamma=0$. In all
examples, the value of $\gamma$ is adaptively changed in every cycle depending
on the smallest negative eigenvalue of (18)

The sparse UKF algorithm is based on the assumption that $P^{a}$ can be
approximated by a sparse matrix $\underline{P}^{a}$. Although the
$\sigma$-points in the algorithm play a similar role as that of ensembles in
EnKF, using sparse UKF one can avoid the problem of rank deficiency. For
systems with very high dimensions, the number of ensemble members used in an
EnKF is much smaller than the dimension. As shown in Figure 1, the narrow and
tall matrix of ensemble vectors makes EnKF fundamentally a rank deficient
approach. In contrast, the block diagonal matrix $\underline{P}^{a}$ shown in
Figure 1 as a sparse approximation of $P^{a}$ has full rank. If $N_{sp}$ of
$\underline{P}^{a}$ is an integer close to the ensemble size of an EnKF, then
the memory required by $\underline{P}^{a}$ is smaller than that by the
ensemble matrix because of the symmetry of error covariance. The required
computational load in Step 1 is extremely high if full state vectors are
computed. Thanks to the sparsity, we only need to compute the entries with
indices in ${\cal I}^{b}$. For a sparse UKF to be successful for problems with
very high dimensions, it is critical to have component-based numerical models
so that subsets of entries defined by ${\cal I}^{b}$ are computed in parallel;
and most entries of state vectors are not evaluated at all. In addition each
component-based computation requires a part of the state vector only. It saves
the computer I/O load.

Figure 1: Patterns of ensemble vectors and sparse error covariances

Memory and (I/O) requirements are big factors influencing the efficiency of
computational algorithms. Because an error covariance and its diagonal blocks
are symmetric, the memory and I/O usage can be reduced by almost a half for
subproblems with symmetry. The granularity of a computational algorithm has
considerable impact on its efficiency in parallel computing. By granularity
control we mean that one can divided a high dimensional problem into
subproblems of desired dimensions. As shown in Figure 1, the error covariance
and its square root consist of sparse blocks or sparse columns. This is
different from the EnKF in which state vectors in an ensemble are dense. In a
sparse approximation of $P^{a}$, the number of nonzero entries, a parameter
similar to that used for distance-based localization methods, can be easily
changed in $\underline{P^{a}}$ so that the error covariance and its square
root can be grouped into smaller blocks of different sizes for parallel
computation.

###  2.3 Lorenz-96 model

In this section, we use a Lorenz-96 model that was first introduced in [6] to
test the performance of the sparse UKF. Consder

$\begin{array}[]{lllllllll}\begin{aligned}
\displaystyle\frac{dx_{i}}{dt}&=(x_{i+1}-x_{i-2})x_{i-1}-x_{i}+F,&i=1,2,\cdots,n,\\\
x_{n+1}&=x_{1},\\\ n&=40,\\\ \Delta t&=0.025,\\\
F&=8.\end{aligned}\end{array}$ (19)

The system has chaotic trajectories as shown in Figure 2, a plot of
$x_{1}(t),x_{2}(t),x_{3}(t)$.

Figure 2: A chaotic trajectory of the Lorenz-96 model, $x_{1}$(solid),
$x_{2}$(dash), $x_{3}$(dot).

The simulations are conducted based on a 4th-order Runge-Kutta discretization.
The trajectories are used as the ground truth. The sparsity pattern for
$\underline{P}^{a}$ and $\sqrt{\underline{P}^{a}}$ are assumed to be centered
along the diagonal line with a fix length. The total number of nonzero entries
in each column is $N_{sp}$. We would like to point out that the sparse
matrices are approximations of the true error covariance and its square root.
The sparsity pattern of $\sqrt{\underline{P}^{a}}$ is, in fact, different from
that of $\underline{P}^{a}$. In the approximation, however, we ignore the
difference and use the same sparsity pattern for both. As a result, the memory
required by $\sqrt{\underline{P}^{a}}$ is reduced almost by half. This idea
works fine for the Lorenz-96 model. However, a systematic way of choosing the
sparsity pattern for $\sqrt{\underline{P}^{a}}$ based on given
$\underline{P}^{a}$ is an open problem that needs further study.

The numerical experimentation is based on $N=1000$ uniformly distributed
random initial states in $[-1\;1]$. The time step size is $\Delta t=0.025$.
The total number of time steps for each simulation is $N_{t}=4000$. The number
of observations at any given time is $m=20$, i.e. every other state variable
is measured. The observational error has the Gaussian distribution with a
covariance $R=I$, the identity matrix. The initial background error covariance
is $P^{b}(0)=0.2I$. The following RSME is used to measure the accuracy of
estimation

$\begin{array}[]{lllllllll}RSME=\sqrt{\displaystyle\frac{1}{n(N_{t}+1)}\displaystyle\sum_{k=0}^{N_{t}}||x^{a}(k)-x^{truth}(k)||^{2}_{2}}.\end{array}$
(20)

For comparison, an EnKF is applied to the same data set. The localization
radius is $\rho=4$ and the inflation factor is $\sqrt{1.08}$. A full scale UKF
based on dense error covariance is also applied as the best estimator for
comparison. As an indication of computational load, the number of entries to
be computed in each algorithm is shown in Table 2. The boxplot of simulation
results is shown in Figure 3. To summarize, the sparse UKF has considerably
smaller error variation than that of the EnKF. This is expected because the
new approach avoids the problem of rank deficiency. The medians of estimation
errors are also smaller than that of EnKF. In the cases of $N_{sp}=7$ and
$11$, the memory size required by the sparse error covariance is smaller than
the memory size needed to store the ensemble vectors in the EnKF if
$N_{ens}=10$. The Cholesky factorization, while maintaining the sparsity
property, may require additional memory. In this example of Lorenz-96, we use
the same sparsity pattern for both $\sqrt{\underline{P}^{a}}$ and
$\underline{P}^{a}$. This assumption simplifies the algorithm and reduces the
memory and I/O requirements. The computational load, in the number of entry
evaluations, is increased for sparse UKF because of the number of
$\sigma$-points is $2n$. Reducing the number of entries being evaluated and
testing the impact of Cholesky factorization on the efficiency of UKF is
ongoing research that is not addressed in this paper.

Filter | Size | Entries | Error | Error | Error
---|---|---|---|---|---
|  | EVAL | Median | Mean | STD
EnKF | $N_{ens}=10$ | 400 | 0.3462 | 1.0741 | 1.0652
S-UKF | $N_{sp}=7$ | 600 | 0.3061 | 0.3067 | 0.0071
S-UKF | $N_{sp}=11$ | 920 | 0.2691 | 0.2691 | 0.0048
Table 2: Summary of simulation results Figure 3: Boxplot of RMSE

##  3 Progressive EKF

In a sparse UKF, the $\sigma$-points are computed by taking a square root of
the error covariance, such as the Cholesky factorization. The resulting
$\sigma$-points are sparse. However, this process may require additional
memory and computation. In this section, we propose a progressive algorithm of
approximating error covariance without taking square roots.

###  3.1 Basic ideas

The main assumption for this algorithm is the following progressive
relationship

$\begin{array}[]{lllllllll}M_{k-1}P^{a}(k-1)M_{k-1}^{T}=P^{a}(k-1)+\Delta
P^{b},\end{array}$ (21)

where $\Delta P^{b}$ is assumed to be small. In (21), $M_{k-1}$ is the
Jacobian of ${\cal M}$ at $x^{a}(k-1)$. Similarly, the Jacobian of ${\cal H}$
is $H_{k}$. To estimate $\Delta P^{b}$, assume

$\begin{array}[]{lllllllll}M_{k-1}&=I+\Delta M_{k-1}.\end{array}$ (22)

where we assume that $\Delta M_{k-1}$ is small. If the system model is based
on the discretization of a differential equation with a small time step size,
then

$\begin{array}[]{lllllllll}{\cal M}(x(k-1))=x(k-1)+O(\Delta
t^{\alpha}),\;\;\alpha\geq 1.\end{array}$ (23)

The Jacobian of $O_{k-1}(\Delta t^{\alpha})$ in space variables is expected to
have small value if $\Delta t$ is small, which makes (22) a reasonable
assumption. Then we have

$\begin{array}[]{lllllllll}\begin{aligned} &M_{k-1}P^{a}(k-1)M_{k-1}^{T}\\\
&=(I+\Delta M_{k-1})P^{a}(k-1)(I+\Delta M_{k-1}^{T})\\\ &=P^{a}(k-1)+\Delta
M_{k-1}P^{a}(k-1)+\left(\Delta M_{k-1}P^{a}(k-1)\right)^{T}+\Delta
M_{k-1}P^{a}(k-1)\Delta M_{k-1}^{T}\\\ &\approx P^{a}(k-1)+\Delta
M_{k-1}P^{a}(k-1)+\left(\Delta
M_{k-1}P^{a}(k-1)\right)^{T}.\end{aligned}\end{array}$ (24)

This is in consistent with (21). It can be computed using a tangent linear
model. Or it can be approximated using the dynamical model

$\begin{array}[]{lllllllll}\begin{aligned} &M_{k-1}P^{a}(k-1)M_{k-1}^{T}\\\
&=(I+\Delta M_{k-1})P^{a}(k-1)(I+\Delta M_{k-1}^{T})\\\ &\approx\left({\cal
M}(x(k-1)+\delta P^{a}(k-1))-{\cal M}(x(k-1))\right)/\delta\\\ &+\left({\cal
M}(x(k-1)+\delta P^{a}(k-1))-{\cal
M}(x(k-1))\right)^{T}/\delta-P^{a}.\end{aligned}\end{array}$ (25)

where $\delta>0$ is the step size of a finite difference approximation of
$\Delta M_{k-1}P^{a}$. Its value should be determined depending on the
numerical model and its linearization. In (25), a vector and matrix summation
is a new matrix resulting from adding the vector to every column in the
matrix. Applying an operator to a matrix is to apply the operator to every
column in the matrix.

###  3.2 Progressive EKF

The column vectors in the matrices in (24) and (25) are sparse. However, the
number of column vectors equals $n$, which can be as high as $10^{6}-10^{7}$
for some atmospheric models. Applying a full model to all the vectors is
impractical because of the high computational and I/O loads. On the other
hand, if we approximate the error covariance using a given sparsity, only a
small portion of the entries in each column vector is to be evaluated.
Evaluating the entire state vector is unnecessary. This is the reason we need
a component-based model. Then the algorithm of progressive EKF is summarized
as follows.

Algorithm II (Progressive EKF)
Given the initial analysis at $t=k-1$,

$\begin{array}[]{lllllllll}x^{a}(k-1)\mbox{ and
}\underline{P}^{a}(k-1).\end{array}$ (26)

Step 1. Forecast

$\begin{array}[]{lllllllll}x^{b}={\cal M}(x^{a}(k-1)),\\\ y^{b}={\cal
H}(x^{b}).\end{array}$ (27)

Step 2. Background error covariance

$\begin{array}[]{lllllllll}\begin{aligned}
\underline{P}^{b}&=\underline{\left({\cal
M}\left(x^{a}(k-1)+\delta\underline{P}^{a}(k-1),{\cal
I}\right)-x^{b}\right)}/\delta\\\ &+\underline{\left({\cal
M}\left(x^{a}(k-1)+\delta\underline{P}^{a}(k-1),{\cal
I}\right)-x^{b}\right)}^{T}/\delta-\underline{P}^{a}+Q.\end{aligned}\end{array}$
(28)

Step 3. Kalman gain and analysis

$\begin{array}[]{lllllllll}\begin{aligned}
K&=\underline{P}^{b}H_{k}^{T}(H_{k}\underline{P}^{b}H^{T}_{k}+R)^{-1},\\\
x^{a}&=x^{b}+K(y_{o}-y^{b}),\\\
\underline{P}^{a}&=(I-KH_{k})\underline{P}^{b}.\\\ \end{aligned}\end{array}$
(29)

Different from the sparse UKF, this algorithm avoids the computation of matrix
square roots. However, the algorithm requires that $\Delta P^{b}$ in (21) can
be approximated effectively. From (23), the method is expected to work better
for a small time step-size. If $\Delta t$ is large, $\Delta M_{k-1}$ in (22)
may not be small enough. A remedy is to use a refined step-size in an inner-
loop computation. More specifically, the discrete model is a discretization of
a continuous-time model. The discrete time moment $k-1$ corresponds to the
continuous time moment $(k-1)\Delta t$. We refine the step size by dividing
the time interval into $n_{p}$ subintervals. In our examples, we choose
$n_{p}=2$. The refined time steps are

$\begin{array}[]{lllllllll}(k-1)\Delta t,(k-1)\Delta
t+\displaystyle\frac{\Delta t}{n_{p}},\cdots,(k-1)\Delta
t+s\displaystyle\frac{\Delta t}{n_{p}},\cdots,k\Delta t,\;\;\;\;\;0\leq s\leq
n_{p}\end{array}$ (30)

For the inner loop, one can compute a sequence of backgrounds,
$\tilde{x}^{b}(s)$.

$\begin{array}[]{lllllllll}t_{s}=(k-1)\Delta t+s\displaystyle\frac{\Delta
t}{n_{p}},\\\ \tilde{x}^{b}(s)=\tilde{\cal
M}_{t_{s}}(x^{a}(k-1)),&s=1,2,\cdots,n_{p}.\end{array}$ (31)

where $\tilde{\cal M}_{t_{s}}$ represents the refined model function in the
time interval from $t=(k-1)\Delta t$ to $t=t_{s}$. In Step 2, repeat (28)
$n_{p}$ times along the sequence of background states, $\tilde{x}^{b}(s)$,
without adding $Q$ until the last round. This refined Step 2 increases the
computational load, while improving the accuracy of the progressive
estimation.

###  3.3 Examples

In the following, we apply the progressive EKF to the Lorenz-96 model using
the same parameters as in (19). The error covariance is approximated using
sparsity matrices with $N_{sp}=7,11,17$. For $N_{sp}=11$, we tested the idea
of refining step-size using $n_{p}=1$ and $n_{p}=2$. The results are shown in
Figure 4 and summarized in Table 3. Comparing to EnKF, the error variations of
the progressive EKF are uniformly and significantly smaller. If $N_{sp}=7$,
which is smaller than the ensemble size $N_{ens}=10$, the median value of
estimation error is larger than that of the EnKF. The median error for
$N_{sp}=11$ is comparable to that of the EnKF. If a refined step-size in (31)
is applied, for instance $n_{p}=2$, the median estimation error is further
reduced. Comparing to the performance of the sparse UKF, the error variations
are similar. However, the estimation error of the sparse UKF has a smaller
median in all cases. For example, to achieve a similar performance as the
sparse UKF when $N_{sp}=11$, one has to use a larger sparsity index
$N_{sp}=17$ for the progressive EKF.

Figure 4: Boxplot of RMSE. For Prograssive-KF, $N_{sp}$ = 7, 11, and 17

The computational load of the progressive EKF, in terms of the component
evaluation, is comparable to the EnKF, as shown in Table 3.

Filter | Size | Entries | Error | Error | Error
---|---|---|---|---|---
|  | EVAL | Median | Mean | STD
EnKF | $N_{ens}=10$ | 400 | 0.3462 | 1.0741 | 1.0652
P-KF | $N_{sp}=7$ |  |  |  |
| $N_{p}=1$ | 320 | 0.3845 | 0.3846 | 0.0055
P-KF | $N_{sp}=11$ |  |  |  |
| $N_{p}=1$ | 480 | 0.3455 | 0.3458 | 0.0050
P-KF | $N_{sp}=11$ |  |  |  |
| $N_{p}=2$ | 480x2 | 0.3041 | 0.3041 | 0.0044
P-KF | $N_{sp}=17$ |  |  |  |
| $N_{p}=3$ | 720x2 | 0.2872 | 0.2873 | 0.0046
Table 3: Summary of simulation results

##  4 Conclusions

Two algorithms of Kalman Filters based on sparse error covariances are
introduced. They are tested using the Lorenz-96 model with $40$ state
variables and chaotic trajectories. Both algorithms share the same basic idea:
the error covariance is approximated using a sparse matrix. Thanks to the
sparsity, the required memory size is significantly reduced. The symmetry of
the error covariance can potentially reduce the I/O load. The analysis error
covariance can be updated as a sparse matrix in each cycle using a
deterministic process, either a square root matrix or a progressive algorithm.
The updated sparse matrix is then used as the background error covariance for
the next cycle. Relative to the EnKF, the main advantage of the proposed
methods is that the estimation process do not need an ensemble; and the error
covariance has a full rank. The algorithms do not suffer issues of rank
deficiency as in EnKFs. As a result, the variation of analysis error is
constantly small in all examples. Techniques of localization and covariance
inflation are unnecessary. Relative to 4D-Var methods, the proposed algorithms
are mostly parallel. They provide not only the state estimate but also the
analysis error covariance. For the purpose of scalability, we suggest that the
proposed methods are applied with component-based numerical models. From the
examples, the sparse UKF has better accuracy than the progressive EKF. If the
computational load of taking square roots of sparse matrices is affordable,
then the sparse UKF is the approach of our choice. On the other hand, the
progressive EKF is a simple algorithm that avoids taking square roots of large
matrices, provided that the progressive approximation of error covariance is
adequately accurate. Although most conclusions drawn in this paper are based
on simulations using the Lorenz-96 model, the algorithms are developed for
general applications. Testing the methods using different types of system
models is a main topic of our future work.

## References

  * [1] T. A. Davis, Direct Methods for Sparse Linear Systems, SIAM, 2006.
  * [2] T. Davis, S. Rajamanickam, and W. M. Sid-Lakhdar, A survey of direct methods for sparse linear systems, Acta Numerica, Vol. 25, 2016, pp. 383-566.
  * [3] P. L. Houtekamer and F. Zhang, Review of the Ensemble Kalman Filter for Atmospheric Data Assimilation, Monthly Weather Review, Vol. 144, 2016, pp. 4489 - 4532.
  * [4] S. Julier, J. Uhlmann, and H. F. Durrant-Whyte, A New Method for the Nonlinear Transformation of Means and Covariances in Filters and Estimators, IEEE Trans. Automatic Control, Vol. 45 (3), 2000, pp. 477-482.
  * [5] S. J. Julier and J. K. Uhlmann, Unscented filtering and nonlinear estimation, Proceedings of the IEEE, Vol. 92, No. 3, 2004, pp.
  * [6] E. Lorenz, Predictability – A problem partly solved. Seminar on Predictability, Vol. I, ECMWF, 1996.
  * [7] E. Rozin and S. Toledo, Locality of reference in sparse Cholesky factorization methods, Electronic Transactions on Numerical Analysis, Vol. 21, 2005, pp. 81-106.
  * [8] L. Xu, R. Rosmond, and R. Daley, Development of NAVDAS-AR: Formulation and initial tests of the linear problem, Tellus, 57A, 2005, pp. 546-559.