# Dense Air Quality Maps Using Regressive Facility Location Based Drive By
Sensing

Charul Paliwal  Indraprastha Institute of Information Technology, Delhi
Pravesh Biyani  Indraprastha Institute of Information Technology, Delhi

###### Abstract

Currently, fixed static sensing is a primary way to monitor environmental data
like air quality in cities. However, to obtain a dense spatial coverage, a
large number of static monitors are required, thereby making it a costly
option. Dense spatiotemporal coverage can be achieved using only a fraction of
static sensors by deploying them on the moving vehicles, known as the drive by
sensing paradigm. The redundancy present in the air quality data can be
exploited by processing the sparsely sampled data to impute the remaining
unobserved data points using the matrix completion techniques. However, the
accuracy of imputation is dependent on the extent to which the moving sensors
capture the inherent structure of the air quality matrix. Therefore, the
challenge is to pick those set of paths (using vehicles) that perform
representative sampling in space and time. Most works in the current
literature for vehicle subset selection focus on maximizing the spatiotemporal
coverage by maximizing the number of samples for different locations and time
stamps which is not an effective representative sampling strategy and results
in noisy imputations. We present regressive facility location based drive by
sensing, an efficient vehicle selection framework that incorporates the
smoothness in neighboring locations and autoregressive time correlation while
selecting the optimal set of vehicles for effective spatiotemporal sampling.
We show that the proposed drive by sensing problem is submodular, thereby
lending itself to a greedy algorithm but with performance guarantees. We
evaluate our framework on selecting a subset from the fleet of public
transport in Delhi, India. We illustrate that the method proposed in the paper
samples the representative spatiotemporal data against the baseline methods,
reducing the extrapolation error on the simulated air quality data. Our
method, therefore, has the potential to provide cost effective dense air
quality maps.

##  1 Introduction

Air pollution has drawn attention in recent years because of its profound
effect on people’s health [1]. It is crucial to assess the pollution level in
a region by monitoring the air quality levels and recommend strategies to
combat pollution. There are significant factors that affect the air quality of
an area, such as transportation, electricity, fuel uses, industrial parameters
such as power plant emissions, etc. Hence, air quality monitoring is required
to make strategies for emission control as well as verifying strategies that
control emissions.

The most widely used method to monitor air quality is by using static sensing,
that is mounting sensors at fixed locations. The temporal resolution of static
sensing is high, i.e., the air quality data is available for almost all the
sampled timestamps. However, the spatial coverage in a region depends on the
number of sensors installed which is generally limited owing to the cost
constraints. We can leverage the spatial and temporal correlation in the air
quality data to perform cost-effective spatiotemporal monitoring. We can
share/multiplex an air quality sensor over multiple locations and perform
sequential temporal sensing without losing air quality information. Moving
sensor paradigm can be a candidate to obtain such dense AQ (Air Quality) map
without needing an expensive setup of hundreds of static monitors [2, 3]. With
the technological progress in AQ sensing, it is now possible to put these
monitors in buses and other public transit and perform spatiotemporal AQ
sensing. However, the moving sensors reading would be spatially and temporally
sparse and incomplete because of the unavailability of the moving vehicle
across all the spatial locations and timestamps.

Thanks to the inherent temporal and spatial redundancy available in the air
quality data, one can potentially impute the remaining (incomplete) data by
utilising the matrix completion techniques [4, 5] and obtain a dense air
quality map in a fraction of the corresponding fixed sensing cost. However,
the quality of imputation crucially relies on the ”quality” of drive by
sensing and the choice of ”routes” that vehicles take to sample the air
quality.

The problem of dense air quality maps using drive by sensing is defined as:
Given a binary bus occupancy data
${\bm{\mathcal{Y}}}\in\\{0,1\\}^{|\mathcal{L}|\times|\mathcal{T}|\times|\mathcal{B}|}$,
where $\mathcal{L}$ denotes the set of locations, $\mathcal{T}$ denotes the
set of timestamps and $\mathcal{B}$ denotes the buses. The entry $y_{l,t,b}=1$
if the bus b samples the location $l$ in the time stamp $t$. The aim is to
pick the best subset of buses $\mathcal{M}\in\mathcal{B}$ that perform
representative sampling for the spatiotemporal AQ matrix. Subsequently, by
appropriately modeling the spatial and temporal structure present in the AQ
data sampled by the set of buses $\mathcal{M}$, one can extrapolate/impute the
missing data in the spatiotemporal matrix. In a nutshell, the aim is to pick
those buses that exploit the spatiotemporal structure in the air quality data
resulting in an effective imputation and thereby creating dense anytime-
anywhere AQ map.

Limited work on selecting the optimal set of vehicles for sensing
spatiotemporal data models each entry in the spatiotemporal matrix as
independent and tries to maximize the number of entries in the matrix [6, 7,
8]. Considering all the locations as independent ignores the correlation
across locations, thereby highlighting the sub-optimality of such an approach.
Indeed, sampling two consecutive neighborhood locations are less effective
than sampling the diverse spread across locations in an area.

Our work exploits the smoothness in the spatial locations to select buses that
sample diverse and representative sets of locations. Further, the air quality
data for a given location is smooth in time, i.e., the sensor data varies
slowly in time with decreasing correlation as the interval increases. Thereby,
sampling the consecutive timestamps of a site is a less effective strategy
than sampling at distant time stamps. The proposed Regressive Facility
Location (RFL) framework encapsulates the temporal smoothness using an
autoregressive time series structure to sample the spatiotemporal data. Also,
since the future temporal data is not available to create the dense map at a
particular timestamp, we use the causal temporal smoothness in the proposed
RFL framework. This ensures that if a location $l$ is sampled at a given time,
other buses can bypass sampling location $l$ and the nearby locations in the
subsequent neighboring future time stamps. The Regressive Facility Location
(RFL) framework proposes to encapsulate the slowly varying temporal data
pattern into the facility location framework and select the buses such that
the sampled locations will be representative of the area and will be diverse
across all time stamps.

Practically it is not possible to access the real-world air quality data
without deploying the sensors on the selected set of buses. Therefore, to
evaluate the performance of the proposed RFL drive by sensing framework, we
simulate real-world air quality data. Further, we obtain the dense air quality
maps using the matrix completion framework from the sampled spatio-temporal
data and observe that RFL provide more accurate dense air quality maps. The
overall contribution of our work are:

  1. 1.

We propose Regressive Facility Location (RFL), a novel method to select a
fleet of vehicles (buses) for drive by sensing that incorporates the
smoothness of the sensor reading across space as well as time to provide
effective data sampling.

  2. 2.

We show that the proposed RFL algorithm gain is submodular, therefore a greedy
heuristic provides a solution within an approximation ratio of
$(1-\frac{1}{e})$.

  3. 3.

We simulate the smooth data over space and time and sample the area using the
proposed drive by sensing framework. To evaluate the performance of the
proposed framework, we create dense maps using a matrix completion framework
from the sampled data.

###  1.1 Notations

A tensor is represented by $\bm{\mathcal{Y}}$. Matrix is denoted by
$\mathbf{Y}$. A set is denoted by $\mathcal{Y}$, the number of element in a
set is denoted by $Y$ where $Y=|\mathcal{Y}|$, the element of the set is
denoted by $y$.

##  2 Background Review

###  2.1 Vehicle Subset Selection

There has been limited work on selecting the optimal set of vehicles for
sensing spatiotemporal data. Authors in [8] propose mobile sensor placement
for vehicles to maximize coverage. Authors in [6] propose a drive by sensing
framework for taxis and buses. Authors in [7] proposed Point of interest
oriented (POIs) Bus Selection algorithm to select buses where the coverage of
a bus is defined in terms of the POIs. Authors in [9] proposed the optimal
placement of reference monitors to make mobile sensors k-hop calibrable.
However, the mobile sensors are fixed and reference monitors are optimally
selected. All these framework maps the problem of vehicle subset selection to
the maximum cover problem or set cover problem to maximize the spatiotemporal
coverage.

###  2.2 Greedy Submodular Maximization

Definition 1: A function f is monotonically non decreasing if
$\forall\,\,\,\mathcal{C}\subseteq\mathcal{D}$

$f(\mathcal{C})\leq f(\mathcal{D})$

Definition 2: A function f: $2^{\mathcal{B}}\rightarrow\mathbf{R}$ is
submodular if
$\forall\,\,\,\mathcal{C}\subseteq\mathcal{D}\subseteq\mathcal{B}$ and every
$b\in\mathcal{B}\backslash\mathcal{D}$

$f(\mathcal{C}\cup b)-f(\mathcal{C})\geq f(\mathcal{D}\cup b)-f(\mathcal{D})$

Given a submodular set function $f$ , maximization of $f$ over all subsets of
size at most $k$ of the ground set $\mathcal{G}$, i.e. $|\mathcal{G}|=k$,

$f(\mathcal{V})=\max\limits_{\mathcal{V}:|\mathcal{V}|\leq k}f(\mathcal{V})$
(1)

is an NP-hard problem [10]. A monotone non decreasing submodular function
solution can be approximate by an greedy algorithm within $(1-\frac{1}{e})$ of
the global maximum [11, 12].

###  2.3 Generating Smooth Graph Signals

To generate a smooth graph signal across nodes, adjacency matrix can be used.
Consider a graph with $\mathcal{L}$ set of nodes/locations, the weighted
adjacency matrix can be defined as $\mathbf{G}\in\mathbf{R}^{L\times L}$ can
be decomposed as.

$\mathbf{G}=\mathbf{U}\Sigma\mathbf{U}^{T}$

The graph Fourier transform of a signal $\bm{x}\in\mathbf{R}^{L}$ is given by
$\hat{\bm{x}}$.

$\hat{\bm{x}}=\mathbf{U}^{-1}\bm{x}$

The actual signal can be expressed as

$\bm{x}=\mathbf{U}\hat{\bm{x}}$

The signal is smooth in the time domain if the signal is bandlimited in the
frequency domain [13, 14, 15, 16, 17] hence,

$\hat{\bm{x}}_{k}=0\,\,\,\forall k\geq m$

Therefore

$\bm{x}=\mathbf{U}_{(m)}\hat{\bm{x}}$ (2)

where $\mathbf{U}_{(m)}\in\mathbf{R}^{L\times m}$ represents the first $m$
eigenvectors of the matrix $\mathbf{G}$.

The signal $\hat{\bm{x}}$ can be sampled from random normal distribution as
described in the literature [14, 17, 16]. Then using the signal $\hat{\bm{x}}$
and $\mathbf{U}$, actual smooth signal can be constructed over the nodes of
the graph via Eq. 2. The signal generated will be smooth on the locations
only. We simulate real-world air quality data that is smooth in space and
time.

###  2.4 Dense maps using Missing data imputation

Air quality data exhibits both spatial and temporal correlation thereby
generating redundancy (not high rank nature) [18, 19]. Low rank matrix
completion has been proposed to estimate the missing spatiotemporal data [4,
20]. Further state space model for incorporating the temporal evolution in the
matrix completion framework is also proposed in the literature [21, 22, 5].
However, dense air quality maps cannot be computed using the traditional
matrix completion framework for the cold start locations. Cold start locations
refers to the set of locations that are never sampled by the selected set of
buses. Similarity matrix can be exploited to impute the data for the cold
start locations [23, 24].

##  3 Bus Selection Framework

Bus selection for drive by sensing is defined as follows: We are given a
binary bus occupancy tensor
${\bm{\mathcal{Y}}}\in\\{0,1\\}^{|\mathcal{L}|\times|\mathcal{T}|\times|\mathcal{B}|}$,
where $\mathcal{L}$ is the set of locations, $\mathcal{T}$ is the set of
timestamps, and $\mathcal{B}$ is the entire set of buses. Buses run on a
specific route with a time schedule. The entry $y_{l,t,b}=1$ denotes the
availability of the sensor reading placed at bus $b$ for a particular location
$l$ and time slot $t$. Our work aims to pick $k$ buses from the entire set of
buses $B$ for drive by sensing, where $\mathcal{M}\subset\mathcal{B}$
represents the selected set of buses. For a subset of buses
$\mathcal{M}\subset\mathcal{B}$, we define a binary sampling matrix
$\bm{\Theta}\in\\{0,1\\}^{L\times T}$ as

$\bm{\Theta}(\mathcal{M})=\max(\sum\limits_{k\in\mathcal{M}}\bm{\mathcal{Y}}_{l,t,k},1)\,\,\,\,\,\forall\,\,(l\in\mathcal{L}),\,\,\,\forall\,\,(t\in\mathcal{T})$
(3)

The distance between two locations $i$ and $j$ are defined as
$\mathbf{D}_{i,j}$. We then define the normalized distance between two
locations $\textbf{ND}_{i,j}=\frac{\mathbf{D}_{i,j}}{\max{(\mathbf{D})}}$. All
the entries in ND are between 0 and 1. The smaller the normalized distance
between two locations, the higher is the similarity. To incorporate the
smoothness in the spatial locations, we define the similarity between
locations in $\mathbf{S}\in\mathbf{R}^{L\times L}$ as:

$\mathbf{S}_{i,j}=1-\textbf{ND}_{i,j}$ (4)

###  3.1 Baseline methods

####  3.1.1 Max Coverage (MC)

We model our problem as the max coverage problem where each entry in the
sampling matrix $\bm{\Theta}$ is considered as an element[8]. We have the set
of buses, each bus $e$ samples some elements in the matrix $\bm{\Theta}(e)$.
Max Coverage framework selects the $k$ number of buses, set
$\mathcal{M}\subset\mathcal{B}$ where $|\mathcal{M}|=k$ such that entries in
the sampling matrix $\bm{\Theta}(\mathcal{M})$ is maximized. The maximum
coverage problem is NP-hard and submodular; there exists a greedy heuristic
that provides a solution within an approximation ratio of
$(1-\frac{1}{e})$[25]. We use Percentage Coverage as a measure for the Max
Coverage Algorithm. The gain for Max Coverage to maximize the occupancy of the
selected buses $\mathcal{M}$ is defined as

$PC(\mathcal{M})=\frac{\sum\limits_{i\in\mathcal{L},j\in\mathcal{T}}\bm{\Theta}(\mathcal{M})_{i,j}}{L\times
T}*100$ (5)

####  3.1.2 Max Coverage over Location (MCL)

We model our problem as the max coverage over location problem where each
location is considered as an element. We have the set of buses, each bus
samples some locations over all the timestamps. Max Coverage over location
framework selects the $k$ number of buses, set $\mathcal{M}\subset\mathcal{B}$
where $|M|=k$ such that maximum locations are sampled. In max coverage over
location, we are maximizing the number of locations sampled by the selected
buses, while in max Coverage framework, we are maximizing the entries of the
sampling matrix $\bm{\Theta}$ to increase the coverage over location as well
as time. We use Percentage Stop coverage as a measure for the Max Coverage
over Location algorithm as shown in Eq. 6.

$PSC(\mathcal{M})=\frac{\sum\limits_{i\in\mathcal{L}}max\\{\sum\limits_{k\in\mathcal{M},j\in\mathcal{T}}\bm{\mathcal{Y}}_{i,j,k},1\\}}{L}*100$
(6)

###  3.2 Facility Location over space (FLS)

Max coverage and max coverage location treat every location as independent
while selecting the set of buses, which sample the spatiotemporal matrix.
However, a correlation exists across the locations readings; the sensor
readings will be smooth over neighboring locations. We thereby use Facility
Location over space (FLS) to model the smoothness in the locations using a
similarity matrix.

For the selected set of buses $\mathcal{M}$, the sampling matrix is defined in
Eq.3. For a time stamp $t$, the sampled locations by the set of buses
$\mathcal{M}$ is defined as $\theta_{t}(\mathcal{M})$,

$\theta_{t}(\mathcal{M})=\\{i\\}\,\,\,\forall\,\,\,\bm{\Theta}_{i,t}(\mathcal{M})|_{=1}$
$FLS(\mathcal{M})=\frac{\sum\limits_{t\in\mathcal{T}}\sum\limits_{l\in\mathcal{L}}\pi^{l}_{t}(\mathcal{M})}{L\times
T}$ (7) $\pi^{l}_{t}(\mathcal{M})=\max\limits_{\forall
m\in\theta_{t}(\mathcal{M})}(\mathbf{S}_{l,m})$ (8)

We compute the maximum pairwise similarity between all the locations
($i\in\mathcal{L}$) and the sampled set locations denoted by
$\theta_{t}(\mathcal{M})$ for a timestamp $t$. The gain defined in Eq. 7 is
maximized when the pairwise similarities between the locations and the nearest
chosen location in the selected set of buses are maximized for all the
timestamps. This ensures that we pick the buses that sample the locations that
are representative of the entire area for all the timestamps.

Theorem 1: The Function defined in Eq. 7 is monotone submodular. Therefore a
greedy heuristic provides a solution within an approximation ratio of
$(1-\frac{1}{e})$.
Proof: The proof is carried out in the appendix.

Figure 1: Spatial and temporal similarity Figure 2: (a) Represents the
locations of an example graph G, (b) The Similarity matrix $S$ for Graph G,
(c)-(d) Sampling matrix $\bm{\Theta}$ for time t=1 and 2, (e)-(f) $\pi_{t}$
for FLS, (g)-(h) $\pi_{t}$ for RFL.

###  3.3 Regressive Facility Location (RFL)

FLS selects the buses that sample representative locations across all
timestamps. However, fails to capture the temporal correlation while selecting
the set of buses and treat all timestamps as independent. The sensor
observations of a location over time are not independent and varies slowly
over time. The FLS framework can be modified to incorporate the space and time
similarity as shown in Fig. 1.

The modified framework incorporates the temporal causal similarity i.e., the
future time data is not available to infer the previous time stamps using
$\rho$. The Facility location incorporating the space and time causal
similarity is called Facility Location over Space-Time (FLST). The similarity
coefficients for location $i$ and timestamp $j$ is denoted in Fig. 1. Blue
edges in the graph denote the similarity between different locations for a
particular timestamp $j$. Note that $\rho=0$ is a special case for FLS, and
only blue edges are used for similarity in the FLS framework. FLST models the
inter-temporal similarity for a particular location using $\rho$ and is
denoted by green edges in Fig 1. The inter location and inter-temporal
similarity are denoted by orange edges. Directed edges represent the
similarity from one node to the other. Note that since we are modeling the
temporal causal similarity, reading at location $i$ is relevant only for the
locations at future timestamps. The gain for the FLST is defined as

$FLST(\mathcal{M})=\frac{\sum\limits_{t\in\mathcal{T}}\sum\limits_{l\in\mathcal{L}}\pi_{t}^{l}(\mathcal{M})}{L\times
T}$ (9)
$\pi_{t}^{l}(\mathcal{M})=\max\limits_{j,m\in\bm{\Theta}(\mathcal{M})}(\mathbf{S}_{l,m}*\mathbf{T}_{t,j})$
(10)

where the Temporal causal similarity is defined as

$\mathbf{T}_{t,j}=\begin{cases}\rho^{t-j},&t\geq j\\\ 0,&t<j\end{cases}$ (11)

It is computationally expensive to compute the Eq. 9 since the comparison for
every location and time is done for all the sampled data $\bm{\Theta}$.
Therefore, we propose a fast algorithm that optimizes the gain defined in Eq.
9, with a reduced computational complexity.

Let $\theta_{t}$ be the set of locations that corresponds to sampled locations
in the sampling matrix $\bm{\Theta}$ for time stamp $t$. For $\pi_{0}=0$,
$\rho\in[0,1]$, the regressive facility location gain is defined by

$RFL(\mathcal{M})=\frac{\sum\limits_{t\in\mathcal{T}}\sum\limits_{l\in\mathcal{L}}\pi_{t}^{l}(\mathcal{M})}{L\times
T}$ (12) $\pi_{t}^{l}(\mathcal{M})=\max((\mathbf{S}_{l,m})_{\forall
m\in\theta_{t}(\mathcal{M})},\rho\pi_{t-1}^{l}(\mathcal{M}))\,\,\,\forall
l\in\mathcal{L}$ (13)

Corollary 1: FLST gain defined via Eq. 28 is equivalent to the RFL gain
defined via Eq. 41.
Proof: The proof is carried out in the appendix.

RFL encapsulates the temporal correlation to calculate the gain for
$\pi_{t}^{l}$ by using $\pi_{t-1}^{l}$. The idea for incorporating the
temporal sampling diversity is if the location is sampled at time $t_{1}$,
then the subsequent neighboring time sampling would not be the best
informative sampling as the data is smooth over time. In RFL, we incorporated
the smoothness over time as shown in Fig. 2. Given that a location $g$ is
sampled at time $t=1$, then the $\pi_{2}$ gain in the next timestamp is lower
in FLS than RFL. Moreover, suppose we have two adjacent locations $a$ and $b$,
that are neighboring nodes with high similarity. In that case, sampling the
location $b$ at timestamp $t_{1}$ will also provide information for the
location $a$ at subsequent neighboring time stamps as shown in Fig. 2. The
greedy approach for the Facility location over space is shown in Algorithm 1.

Theorem 2: The Function defined in Eq. 40 is monotone submodular. Therefore a
greedy heuristic provides a solution within an approximation ratio of
$(1-\frac{1}{e})$.
Proof: The proof is carried out in the appendix.

Computational Complexity: For each $l\in\mathcal{L}$ and $t\in\mathcal{T}$
computing $\pi_{t}^{l}(\mathcal{M})$ via Eq. 8 requires $\Lambda$ comparisons
where $\Lambda<L$ hence the computational complexity for each iteration of the
greedy algorithm for all the buses in FLS is $O(\Lambda LTB)$. For each
$l\in\mathcal{L}$ and $t\in\mathcal{T}$ computing $\pi_{t}^{l}(\mathcal{M})$
via Eq. 28 requires $\Lambda T$ comparisons, hence a cost of $O(\Lambda
LT^{2}B)$ for FLST. We reduce the computation complexity of RFL by computing
$\pi_{t}^{l}(\mathcal{M})$ via Eq. 41 to $O(\Lambda LTB)$.

Algorithm 1 Regressive Facility Location

Input: $k,\,\,\,\mathcal{B},\,\,\,\bm{\mathcal{Y}},\,\,\rho,\,\,\,\mathbf{S}$
Output: $\mathcal{M}\subset\mathcal{B}$ of size $k$

Initialization:$\mathcal{M}\leftarrow\phi$; $\pi_{0}=0$

1: for $i$ =1 to $k$ do

2: for each $e\in\mathcal{B}$ do

3: $\bm{\Theta}(\mathcal{M}\cup e)=\max(\sum\limits_{k\in\mathcal{M}\cup
e}\bm{\mathcal{Y}}_{i,j,k},1)$

4: for each $t\in\mathcal{T}$ do

5: $\theta_{t}(\mathcal{M}\cup
e)=\\{i\\}\,\,\,\forall\,\,\,\bm{\Theta}_{i,t}(\mathcal{M}\cup e)|_{=1}$

6: Compute $\pi_{t}^{l}(\mathcal{M}\cup e)$ using Eq. 41 $\,\,\,(\forall
l\in\mathcal{L})$

7: end for

8: $f(\mathcal{M}\cup e)=\frac{1}{L\times
T}\sum_{t=1}^{T}\sum_{l=1}^{L}\pi_{t}^{l}(\mathcal{M}\cup e)$

9: end for

10: $e^{*}=\arg\max_{e}f(\mathcal{M}\cup e)$

11: $\mathcal{M}\leftarrow\mathcal{M}\cup\\{e^{*}\\}$

12: end for

##  4 Experimentation

###  4.1 Simulating Real World AQ data

####  4.1.1 Similarity matrix

To simulate the real world AQ data on the set of location $\mathcal{L}$ and
timestamps $\mathcal{T}$, we learn the similarity matrix based on the static
real world of Delhi [26]. It contains the AQ data for 33 locations for 3
months. We define the entry of the Similarity matrix
$\mathbf{G}_{i,j}=\exp^{(-\lambda d_{i,j})}$ where $d_{i,j}$ is the distance
between the two locations $i$ and $j$. We learn the parameter $\lambda$ using
linear regression and it is observed to be 0.07676 for the AQ data [26]. We
then compute the similarity matrix $\mathbf{G}$ for the set of location
$\mathcal{L}$ based on the actual distance between locations $d_{i,j}$. We
learn the temporal similarity matrix from the data [26] and use it as a
temporal similarity matrix $\mathbf{H}$.

(a) RFL vs MC for atleast 1 time stamp

(b) RFL vs MC for atleast 10 time stamp

(c) RFL vs MCL for atleast 1 time stamp

(d) RFL vs MCL for atleast 10 time stamp

(e) RFL vs FLS for atleast 1 time stamp

(f) RFL vs FLS for atleast 10 time stamp

Figure 3: Coverage plots for all the locations of Delhi for selected bus
$k$=30. Red points denote the point not sampled by both RFL and comparison
Framework, Pink points denote the points sampled by both RFL and comparison
Framework, Blue points denote the points sampled by RFL but not the comparison
framework, Green points denote the points sampled by comparison framework but
not the RFL

To generate a spatiotemporal matrix $\mathbf{Y}$ that varies smoothly over
space and time, we use two variants described in the following sections.

####  4.1.2 Simulated Data 1

  * •

We first factorize the actual matrix $\mathbf{Z}\in\mathbf{R}^{L\times T}$as

$\mathbf{Z}=\mathbf{A}\mathbf{B}^{T}$

where $\mathbf{A}\in\mathbf{R}^{L\times r}$ and
$\mathbf{B}\in\mathbf{R}^{T\times r}$.

  * •

We construct the matrix $\mathbf{A}$ for the locations using the $\mathbf{G}$
and a bandlimited signal. The similarity matrix $\mathbf{G}$ can be decomposed
as

$\mathbf{G}=\mathbf{U}\Sigma\mathbf{U}^{T}$

We generate the matrix $\mathbf{A}$ with rank $r$ as defined in Eq. 2 as.

$\mathbf{a}_{i}=\mathbf{U}_{(m)}\hat{\mathbf{a}_{i}}\,\,\,\,(\forall
i=1\,\,\,\text{to}\,\,\,r)$

where $a_{i}$ represents the $i^{th}$ column of matrix $\mathbf{A}$, $a_{i}$
is sampled from a random normal distribution with standard deviation (0.5) .

  * •

Similarly we construct the matrix $\mathbf{B}$ for time stamps using
$\mathbf{H}$ and a bandlimited signal.

$\mathbf{H}=\mathbf{V}\Lambda\mathbf{V}^{T}$
$\mathbf{b}_{i}=\mathbf{V}_{(n)}\hat{\mathbf{b}_{i}}\,\,\,\,(\forall
i=1\,\,\,\text{to}\,\,\,r)$

  * •

The overall data $\mathbf{Y}$ is generated using $\mathbf{A}$, $\mathbf{B}$
and noise signal $n$ distributed as zero mean and std of 0.001.

$\mathbf{Y}=\mathbf{A}\mathbf{B}^{T}+\mathbf{N}$

We run the experiments for randomly generated 30 spatiotemporal matrix
$\mathbf{Y}$, where $m$ is randomly chosen from 5 to 15, $n$ is chosen from 5
to 15, $r$ is chosen randomly from 20 to 30.

####  4.1.3 Simulated Data 2

To generate a smooth signal over space and time we used the framework
described in paper [27].

$\mathbf{a}_{t}=\mathbf{U}_{(m)}\hat{\mathbf{a}_{t}}$ (14)

The observed signal can be generated by the following Eqs:

$\mathbf{y}_{t}=\mathbf{z}_{t}+\mathbf{n}_{t}$ (15)
$\mathbf{z}_{t}=\mathbf{R}\mathbf{z}_{t-1}+\mathbf{a}_{t}$ (16)

The state transition matrix $\mathbf{R}$ is defined as a general diagonal
matrix $\mathbf{R}=diag(c_{1},c_{2},\dots c_{T})$, where each $c$ represents
the autocorrelation coefficient that describes the time correlation of the
data with the delayed (one time lag) data.We randomly generate 30
spatiotemporal matrix with $c$ as 1, $m$ is randomly chosen from 5 to 15 and
$r$ is chosen randomly from 20 to 30.

| $k=20$ |  $k=35$ | $k=50$
---|---|---|---
| PSC | PC | FLS | RFL | PSC | PC | FLS | RFL | PSC | PC | FLS | RFL
Random bus | 56.796 | 4.612 | 90.202 | 91.936 | 59.102 | 6.593 | 90.379 | 91.563 | 64.078 | 10.584 | 92.515 | 93.579
Max Coverage | 55.218 | 7.357 | 91.867 | 93.163 | 68.932 | 11.983 | 93.801 | 94.693 | 79.49 | 16.117 | 95.089 | 95.764
Max Cov Loc | 80.704 | 5.141 | 90.719 | 92.593 | 92.233 | 7.756 | 92.734 | 94.143 | 96.359 | 10.293 | 93.935 | 94.941
FLS | 60.316 | 6.279 | 93.164 | 94.159 | 73.908 | 10.685 | 94.924 | 95.673 | 75.728 | 14.482 | 95.804 | 96.363
RFL($\rho=0.95$) | 58.981 | 6.349 | 93.122 | 94.201 | 71.845 | 10.628 | 94.885 | 95.654 | 75.85 | 14.629 | 95.785 | 96.392
RFL($\rho=0.98$) | 60.68 | 6.504 | 93.16 | 94.324 | 72.816 | 10.761 | 94.875 | 95.724 | 79.369 | 14.759 | 95.754 | 96.414
RFL($\rho=0.99$) | 64.442 | 6.394 | 92.842 | 94.223 | 73.908 | 10.68 | 94.812 | 95.691 | 79.126 | 14.525 | 95.665 | 96.385
RFL($\rho=1$) | 71.845 | 4.949 | 91.019 | 92.795 | 86.529 | 8.628 | 93.375 | 94.614 | 93.204 | 12.141 | 94.472 | 95.443
TABLE I: Performance Comparison for selecting $k$ buses | Simulated Data 1 | Simulated Data 2
---|---|---
$k$ | Random | MC | MCL | FLS | RFL95 | RFL98 | RFL1 | Random | MC | MCL | FLS | RFL95 | RFL98 | RFL1
20 | 70.595 | 39.333 | 55.582 | 31.71 | 31.395 | 28.507 | 55.577 | 70.239 | 40.13 | 40.225 | 27.68 | 28.03 | 25.85 | 43.70
35 | 47.062 | 22.29 | 27.938 | 10.69 | 11.065 | 10.664 | 21.209 | 40.26 | 19.906 | 16.293 | 8.395 | 9.433 | 8.563 | 13.75
50 | 31.6 | 14.193 | 16.62 | 7.465 | 7.367 | 6.983 | 13.997 | 26.674 | 11.699 | 10.822 | 5.228 | 5.62 | 4.251 | 8.126
75 | 23.546 | 11.275 | 11.296 | 5.312 | 5.07 | 5.065 | 10.832 | 19.46 | 8.448 | 7.197 | 2.736 | 2.627 | 2.532 | 7.263
100 | 18.844 | 10.375 | 8.65 | 4.289 | 4.271 | 4.112 | 7.517 | 14.732 | 8.035 | 5.785 | 1.973 | 1.817 | 1.79 | 4.487
TABLE II: MRE for Dense Map using VBMC(CS)  | Simulated Data 1 | Simulated Data 2
---|---|---
$k$ | Random | MC | MCL | FLS | RFL95 | RFL98 | RFL1 | Random | MC | MCL | FLS | RFL95 | RFL98 | RFL1
20 | 90.788 | 37.538 | 73.922 | 36.585 | 36.749 | 33.602 | 76.477 | 77.619 | 35.275 | 53.419 | 30.424 | 31.479 | 28.864 | 57.799
35 | 44.446 | 21.008 | 26.995 | 10.934 | 11.609 | 10.124 | 19.386 | 39.429 | 20.8 | 19.972 | 9.817 | 11.193 | 9.297 | 14.858
50 | 27.766 | 13.303 | 15.125 | 7.018 | 7.492 | 6.574 | 11.678 | 28.65 | 13.526 | 11.627 | 5.801 | 6.398 | 4.904 | 8.027
75 | 20.34 | 10.068 | 10.15 | 4.528 | 4.586 | 4.52 | 8.842 | 21.873 | 9.704 | 7.38 | 2.581 | 2.591 | 2.581 | 5.178
100 | 16.33 | 9.52 | 7.187 | 3.809 | 3.769 | 3.684 | 6.385 | 16.922 | 9.201 | 5.243 | 1.76 | 1.638 | 1.524 | 3.547
TABLE III: MRE for Dense Map using VBSF(CS)

###  4.2 Creating Dense Maps from Sampled Data

Dense AQ map can be created by first selecting the $\mathcal{M}$ set of buses
then sampling the data $\mathbf{\hat{Y}}$ from the simulated ground truth AQ
matrix $\mathbf{Y}$ based on the sampling matrix $\bm{\Theta}(\mathcal{M})$.
The final step is to impute the missing data in the matrix $\mathbf{\hat{Y}}$.
To evaluate the performance, we compute the MRE score for the missing data
imputation. There exist a problem of cold start while imputation, i.e., some
of the locations are not sampled at all by the selected set of buses.
Therefore, we use an extended version of the matrix imputation method that
handles cold start cases. The matrix completion framework (VBMC) proposed in
paper [20] imposes a low rank structure on the data to impute the missing data
as:

$\displaystyle\mathcal{L}_{1}=\min_{\mathbf{A},\mathbf{B}}||\mathbf{P}_{\Omega}(\mathbf{Y}-\mathbf{A}\mathbf{B}^{T})||_{{F}}$
(17)

where $\mathbf{A}\in\mathcal{R}^{L\times r}$ and
$\mathbf{B}\in\mathcal{R}^{T\times r}$ and $r=$rank$(\mathbf{Y})<<min(L,T)$
Further, VBSF [5] add a regularization on the matrix $\mathbf{B}$ to
incorporate the temporal evolution in addition to the low rankness (Eq. 52) as

$\displaystyle\mathcal{R}(\mathbf{B})=\sum_{i=1}^{T}||\mathbf{b}_{i}-\mathbf{F}\mathbf{b}_{i-1}||$
(18)

However, these framework does not incorporate for the cold start locations.
Incorporating the similarity matrix $\mathbf{G}$ along with the low rank
matrix completion framework can tackle the cold start problem [23, 24].

$\displaystyle\mathcal{L}_{2}=\min_{\mathbf{A},\mathcal{C}}||(\mathbf{G}-\mathbf{A}\mathcal{C}^{T})||_{{F}}$
(19)

We use Variational Bayesian Matrix Completion (Cold start) and Variational
Bayesian Subspace Filtering (Cold start) as the extended matrix completion
frameworks that handles the cold start location data imputation to evaluate
the performance of the dense AQ maps. We impute the missing data using
VBMC(CS) where we optimize the Eqs. (52, 54). We also impute the missing data
using VBSF(CS) where we optimize the Eqs. (52, 53, 54).

###  4.3 Dataset

We use the Delhi Bus GTFS data [28] to obtain the bus occupancy tensor defined
by ${\bm{\mathcal{Y}}}\in\\{0,1\\}^{L\times T\times B}$

  * •

Locations ($L$): There are total 3210 bus stops. We sample $L$ stops such that
the minimum distance between the stops should be $d$. The total number of
locations $L$ is 824 when $d=500$ m.

  * •

Time stamps (T): We use 10 min sampling from 6 am to 10 pm, resulting in
$T$=96 as the buses run during this time.

  * •

Buses ($B$): The number of buses $B$ is 1476.

The entries in $y_{l,t,b}=1$ if bus $b$ is in the 500 meter radius of location
$l$ for the time stamp $t$.

###  4.4 Evaluation Metrics

The evaluation metrics used for the performance comparison are defined as

  * •

Percentage Stop coverage(PSC) is defined in Eq. (6).

  * •

Percentage Coverage (PC) is defined in Eq. (5)

  * •

FLS gain defined in Eq. 7

  * •

RFL gain defined in Eq. 40 for $\rho=0.98$

  * •

Mean Relative Error (MRE):
$\frac{\left\|\mathbf{Y}-\hat{\mathbf{Y}}\right\|}{\left\|\mathbf{Y}\right\|}$,
where $\hat{\mathbf{Y}}$ is the estimate of $\mathbf{Y}$.

Mean Relative Error is used to evaluate the performance of imputed dense AQ
maps.

###  4.5 Performance Comparison

####  4.5.1 Effect of $\rho$ in RFL:

As the $\rho$ increases, PSC increases. For $\rho=1$, the performance is
comparable to the MCL, since if a location is sampled for a timestamp, the
gain corresponding to sampling other timestamps will be zero; therefore MCL
will select the buses that sample different locations to increase the
Percentage stop coverage (PSC). However, the performance of RFL ($\rho=1$)
worsens for PC, FLS and RFL. For $\rho=0.95$, the performance is comparable to
the FLS as the effect of temporal correlation decreases with time in RFL for
lower values of $\rho$. Since the FLS does not incorporate the temporal
correlation in the framework and can sample the same set of representative
locations across timestamps, we observe a reduced PSC for RFL ($\rho=0.95$)
and FLS. We observe similar performance for RFL($\rho=0.98$) and
RFL($\rho=0.99$). However we observe an improved FLS score for
RFL($\rho=0.98$) through out as compared to RFL($\rho=0.99$). Therefore we
evaluate the performance of RFL for $\rho=0.98$ rather than $\rho=0.99$ for
dense map creation.

####  4.5.2 Coverage Plot:

To demonstrate the representations in space and time, we plot the coverage
plot as shown in Fig 3. Every point on the plot is one of the bus stop
locations. Fig 3 (a-c) denotes the locations that are sampled for atleast one
timestamp, thereby illustrating the overall stops coverage. Fig 3 (d-f)
represents the locations that are sampled for atleast 10 timestamps, thereby
demonstrating the temporal coverage for the sampled stops. Pink points denote
the locations sampled by both the RFL and the compared framework. Blue points
denote the point sampled by the RFL but not the compared framework, green
points denote the points sampled by the comparison framework but not the RFL
and red points denote the points not sampled by RFL and the compared
framework. From Fig. 3(a,d), it is observed that the MC does not sample a
diverse set of locations. MCL provides improved spatial coverage, however the
temporal coverage is worse as compared to RFL as shown in Fig. 3(b,e). It can
also be observed that RFL samples a diverse set of locations than FLS as shown
in Fig 3(c,f).

####  4.5.3 Dense AQ map MRE:

The MRE scores for creating dense AQ map is shown in Table II, III. For
simplicity we denote the RFL($\rho=0.95$) as RFL95. It is observed that
RFL($\rho=0.98$) outperforms all the baseline algorithms for dense map
creation VBMC(CS) and VBSF(CS). Since the generated data have temporal
correlation thereby VBSF(CS) performance is better than VBMC(CS) for most of
the cases, especially for higher bus sampling. However, for the lower sampling
of buses (20), it is observed that VBMC(CS) performs better as there are
significantly lower samples to learn the temporal pattern along with the low
rankness.

All the experiments are run on Matlab/Python with the system configuration of
2.3 GHz and 16 GB RAM.

##  5 Conclusion

The anytime-anywhere AQ map is the holy grail of AQ monitoring and is crucial
for combating air pollution, especially in the developing world. To obtain a
dense air quality map, the cities may require only a small fraction of moving
sensors as compared to all static sensors setup. The spatiotemporal
correlation in the air quality data can be leveraged to select the moving
sensor for sampling and to facilitate an effective spatiotemporal
extrapolation on the sampled data. This paper proposes a Regressive Facility
Location for drive-by sensing to select the set of buses that samples the
spatiotemporal AQ data. It is shown that the chosen set of buses provides
representative coverage across space and time. We further obtain the dense AQ
maps using the matrix completion framework from the sampled spatiotemporal
data and observe that RFL provides more accurate dense AQ maps.

##  6 Acknowledgement

The authors would like to thank Dr. Rishabh Iyer (University of Texas, Dallas)
for his suggestions.

## References

  * [1] P. Florina and H. W. Ven-dee, “The global distribution of air pollution,” https://datatopics.worldbank.org/world-development-indicators/stories/the-global-distribution-of-air-pollution.html, 2019\.
  * [2] A. Anjomshoaa, F. Duarte, D. Rennings, T. J. Matarazzo, P. deSouza, and C. Ratti, “City scanner: Building and scheduling a mobile sensing platform for smart city services,” _IEEE Internet of things Journal_ , vol. 5, no. 6, pp. 4567–4579, 2018.
  * [3] D. Hasenfratz, O. Saukh, C. Walser, C. Hueglin, M. Fierz, T. Arn, J. Beutel, and L. Thiele, “Deriving high-resolution urban air pollution maps using mobile sensor nodes,” _Pervasive and Mobile Computing_ , vol. 16, pp. 268–285, 2015.
  * [4] M. T. Asif, N. Mitrovic, J. Dauwels, and P. Jaillet, “Matrix and tensor based methods for missing data estimation in large traffic networks,” _IEEE Trans. Intell. Transp. Syst._ , vol. 17, no. 7, pp. 1816–1825, 2016.
  * [5] C. Paliwal, U. Bhatt, P. Biyani, and K. Rajawat, “Traffic estimation and prediction via online variational bayesian subspace filtering,” _IEEE Transactions on Intelligent Transportation Systems_ , 2021.
  * [6] D. Agarwal, S. Iyengar, M. Swaminathan, E. Sharma, A. Raj, and A. Hatwar, “Modulo: Drive-by sensing at city-scale on the cheap,” in _Proceedings of the 3rd ACM SIGCAS Conference on Computing and Sustainable Societies_ , 2020, pp. 187–197.
  * [7] Y. Gao, W. Dong, K. Guo, X. Liu, Y. Chen, X. Liu, J. Bu, and C. Chen, “Mosaic: A low-cost mobile sensing system for urban air quality monitoring,” in _IEEE INFOCOM 2016-The 35th Annual IEEE International Conference on Computer Communications_. IEEE, 2016, pp. 1–9.
  * [8] J. Ali and V. Dyo, “Coverage and mobile sensor placement for vehicles on predetermined routes: A greedy heuristic approach,” 2017.
  * [9] K. Fu, W. Ren, and W. Dong, “Multihop calibration for mobile sensing: K-hop calibratability and reference sensor deployment,” in _IEEE INFOCOM 2017-IEEE Conference on Computer Communications_. IEEE, 2017, pp. 1–9.
  * [10] S. Fujishige, _Submodular functions and optimization_. Elsevier, 2005.
  * [11] G. Calinescu, C. Chekuri, M. Pál, and J. Vondrák, “Maximizing a submodular set function subject to a matroid constraint,” in _International Conference on Integer Programming and Combinatorial Optimization_. Springer, 2007, pp. 182–196.
  * [12] M. Feldman, J. Naor, and R. Schwartz, “A unified continuous greedy algorithm for submodular maximization,” in _2011 IEEE 52nd Annual Symposium on Foundations of Computer Science_. IEEE, 2011, pp. 570–579.
  * [13] A. Anis, A. Gadde, and A. Ortega, “Efficient sampling set selection for bandlimited graph signals using graph spectral proxies,” _IEEE Transactions on Signal Processing_ , vol. 64, no. 14, pp. 3775–3789, 2016.
  * [14] S. Chen, R. Varma, A. Sandryhaila, and J. Kovavcevi’c, “Discrete signal processing on graphs: Sampling theory¡? pub _newline=””?” _IEEE transactions on signal processing_ , vol. 63, no. 24, pp. 6510–6523, 2015.
  * [15] S. Chen, R. Varma, A. Singh, and J. Kovavcevi’c, “Signal recovery on graphs: Fundamental limits of sampling strategies,” _IEEE Transactions on Signal and Information Processing over Networks_ , vol. 2, no. 4, pp. 539–554, 2016.
  * [16] A. Hashemi, R. Shafipour, H. Vikalo, and G. Mateos, “Accelerated sampling of bandlimited graph signals,” _arXiv preprint arXiv:1807.07222_ , 2018.
  * [17] A. Jayawant and A. Ortega, “A distance-based formulation for sampling signals on graphs,” in _2018 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_. IEEE, 2018, pp. 6318–6322.
  * [18] P. D. Sampson, A. A. Szpiro, L. Sheppard, J. Lindström, and J. D. Kaufman, “Pragmatic estimation of a spatio-temporal air quality model with irregular monitoring data,” _Atmospheric Environment_ , vol. 45, no. 36, pp. 6593–6606, 2011.
  * [19] W. Min and L. Wynter, “Real-time road traffic prediction with spatio-temporal correlations,” _Transportation Research Part C: Emerging Technologies_ , vol. 19, no. 4, pp. 606–616, 2011.
  * [20] S. D. Babacan, M. Luessi, R. Molina, and A. K. Katsaggelos, “Sparse bayesian methods for low-rank matrix estimation,” _IEEE Transactions on Signal Processing_ , vol. 60, no. 8, pp. 3964–3977, 2012.
  * [21] J. Luttinen, “Fast variational Bayesian linear state-space model,” in _Joint European Conference on Machine Learning and Knowledge Discovery in Databases_. Springer, 2013, pp. 305–320.
  * [22] H.-F. Yu, N. Rao, and I. S. Dhillon, “Temporal regularized matrix factorization for high-dimensional time series prediction,” in _Advances in neural information processing systems_ , 2016, pp. 847–855.
  * [23] T. Zhou, H. Shan, A. Banerjee, and G. Sapiro, “Kernelized probabilistic matrix factorization: Exploiting graphs and side information,” in _Proceedings of the 2012 SIAM international Conference on Data mining_. SIAM, 2012, pp. 403–414.
  * [24] I. Barjasteh, R. Forsati, D. Ross, A.-H. Esfahanian, and H. Radha, “Cold-start recommendation with provable guarantees: A decoupled approach,” _IEEE Transactions on Knowledge and Data Engineering_ , vol. 28, no. 6, pp. 1462–1474, 2016.
  * [25] D. S. Hochbaum, “Approximating covering and packing problems: set cover, vertex cover, independent set, and related problems,” in _Approximation algorithms for NP-hard problems_ , 1996, pp. 94–143.
  * [26] CPCB, “Delhi real time air quality data,” 2019. [Online]. Available: https://cpcb.nic.in/real-time-air-quality-data
  * [27] Y. Liu, W. Guo, K. You, L. Zhao, T. Peng, and W. Wang, “Graph learning for spatiotemporal signals with long-and short-term characterization,” _IEEE Transactions on Signal and Information Processing over Networks_ , vol. 6, pp. 699–713, 2020.
  * [28] OTD-Delhi, “Open transit data - delhi,” 2020. [Online]. Available: https://opendata.iiitd.edu.in

##  7 Appendix

###  7.1 Proof of Theorem 1

Let $\theta_{t}$ be the set of locations that corresponds sampled locations in
the sampling matrix $\bm{\Theta}$ for time stamp (column) $t$. The FLS gain is
defined as:

$FLS(\mathcal{M})=\frac{\sum\limits_{t\in\mathcal{T}}\sum\limits_{l\in\mathcal{L}}\pi^{l}_{t}(\mathcal{M})}{L\times
T}$ (20)

where

$\pi^{l}_{t}(\mathcal{M})=\max\limits_{\forall
m\in\theta_{t}(\mathcal{M})}(\mathbf{S}_{l,m})$ (21)

####  7.1.1 Monotone non decreasing

A function f is submodular if $\forall\,\,\,\mathcal{C}\subseteq\mathcal{D}$

$f(\mathcal{C})\leq f(\mathcal{D})$

Proof:
Let $\mathcal{C}$ and $\mathcal{D}$ be the subset of buses following
$\forall\,\,\,\mathcal{C}\subseteq\mathcal{D}\subseteq\mathcal{B}$. We define
$\bm{\Theta}(\mathcal{C})$ and $\bm{\Theta}(\mathcal{D})$ as the sampling
matrix for the subset of buses $\mathcal{C}$ and $\mathcal{D}$.
Let $\theta_{t}$ be the set of locations that corresponds to 1 in the sampling
matrix $\bm{\Theta}(\mathcal{C})$ for the time stamp $t$ and Let
$\theta_{t}^{{}^{\prime}}$ be the set of locations that corresponds to 1 in
the sampling matrix $\bm{\Theta}(\mathcal{D})$ for the time stamp $t$.

$\theta_{t}=\\{i\\}\,\,\,\forall\,\,\,\bm{\Theta}_{i,t}(\mathcal{C})|_{=1}$
$\theta_{t}^{{}^{\prime}}=\\{i\\}\,\,\,\forall\,\,\,\bm{\Theta}_{i,t}(\mathcal{D})|_{=1}$

Since $\mathcal{C}\subseteq\mathcal{D}$, therefore
$\theta_{t}\subseteq\theta_{t}^{{}^{\prime}}$
$\forall\,\,\,t=1\,\,\,\text{to}\,\,\,T$.

For all $l\in\mathcal{L}$, $t\in\mathcal{T}$ we have

$\max\limits_{\forall
m\in\theta_{t}}(\mathbf{S}_{l,m})\leq\max\limits_{\forall
m\in\theta_{t}^{{}^{\prime}}}(\mathbf{S}_{l,m})$ (22)

Hence $f(\mathcal{C})\leq f(\mathcal{D})$

####  7.1.2 Submodular

A function f is submodular if
$\forall\,\,\,\mathcal{C}\subseteq\mathcal{D}\subseteq\mathcal{B}$ and
$b\in\mathcal{B}\backslash\mathcal{D}$

$f(\mathcal{C}\cup b)-f(\mathcal{C})\geq f(\mathcal{D}\cup b)-f(\mathcal{D})$

Consider a bus $b\in\mathcal{B}\backslash\mathcal{D}$, We add the bus $b$ in
both $\mathcal{C}$ and $\mathcal{D}$ and analyse the gain. Let the bus
traverse the set of location $G_{t}$ for a time stamp $t$.

$G_{t}=\\{i\\}\,\,\,\forall\,\,\,\bm{\Theta}_{i,t}(b)|_{=1}$

For all $l\in L$ and $t=1\,\,\,\text{to}\,\,\,T$ we have

$\displaystyle\pi_{t}^{l}(\mathcal{C}\cup
b)-\pi_{t}^{l}(\mathcal{C})=\max\limits_{\forall m\in\theta_{t}\cup
G_{t}}(\mathbf{S}_{l,m})-\max\limits_{\forall
m\in\theta_{t}}(\mathbf{S}_{l,m})$
$\displaystyle=\max(0,\,\,\,\max\limits_{\forall m\in
G_{t}}(\mathbf{S}_{l,m})-\max\limits_{\forall
m\in\theta_{t}}(\mathbf{S}_{l,m}))$ (23)

Since

$\max\limits_{\forall
m\in\theta_{t}}(\mathbf{S}_{l,m})\leq\max\limits_{\forall
m\in\theta_{t}^{{}^{\prime}}}(\mathbf{S}_{l,m})$

Therefore,

$\displaystyle\pi_{t}^{l}(\mathcal{C}\cup
b)-\pi_{t}^{l}(\mathcal{C})=\max(0,\,\,\,\max\limits_{\forall m\in
G_{t}}(\mathbf{S}_{l,m})-\max\limits_{\forall
m\in\theta_{t}}(\mathbf{S}_{l,m}))$
$\displaystyle\geq\max(0,\,\,\,\max\limits_{\forall m\in
G_{t}}(\mathbf{S}_{l,m})-\max\limits_{\forall
m\in\theta_{t}^{{}^{\prime}}}(\mathbf{S}_{l,m}))$ (24)

$\pi_{t}^{l}(\mathcal{C}\cup
b)-\pi_{t}^{l}(\mathcal{C})\geq\max(0,\,\,\,\max\limits_{\forall m\in
G_{t}}(\mathbf{S}_{l,m})-\max\limits_{\forall
m\in\theta_{t}^{{}^{\prime}}}(\mathbf{S}_{l,m}))$ (25)
$\pi_{t}^{l}(\mathcal{C}\cup
b)-\pi_{t}^{l}(\mathcal{C})\geq\max\limits_{\forall
m\in\theta_{t}^{{}^{\prime}}\cup G_{t}}(\mathbf{S}_{l,m})-\max\limits_{\forall
m\in\theta_{t}^{{}^{\prime}}}(\mathbf{S}_{l,m})$ (26)
$\pi_{t}^{l}(\mathcal{C}\cup
b)-\pi_{t}^{l}(\mathcal{C})\geq\pi_{t}^{l}(\mathcal{D}\cup
b)-\pi_{t}^{l}(\mathcal{D})$ (27)

From Eqs (23-27) for all $l\in\mathcal{L}$ and
$t=1\,\,\,\text{to}\,\,\,\mathcal{T}$ we have,

$f(\mathcal{C}\cup b)-f(\mathcal{C})\geq f(\mathcal{D}\cup b)-f(\mathcal{D})$

###  7.2  Prove for Corollary 1

FLST gain defined via Eq. 28 is equivalent to the RFL gain defined via Eq. 41.

$\lambda_{t}^{l}(\mathcal{M})=\max\limits_{m,j\in\bm{\Theta}(\mathcal{M})}(\mathbf{S}_{l,m}*\mathbf{T}_{t,j})$
(28) $\mathbf{T}_{t,j}=\begin{cases}\rho^{t-j},&t\geq j\\\ 0,&t<j\end{cases}$
(29) $\pi_{t}^{l}(\mathcal{M})=\max((\mathbf{S}_{l,m})_{\forall
m\in\theta_{t}(\mathcal{M})},\rho\pi_{t-1}^{l}(\mathcal{M}))\,\,\,\forall
l\in\mathcal{L}$ (30)

Proof:
For all $l\in\mathcal{L}$ and $t=1$ using Eq. 28 we have,

$\lambda_{1}^{l}(\mathcal{M})=\max(\max(\mathbf{S}_{l,m})_{\forall
m\in\theta_{1}(\mathcal{M})},V)\,\,\,\forall l\in\mathcal{L}$ (31)
$V=\max(\mathbf{S}_{l,m})_{\forall
m\in\theta_{j}(\mathcal{M})}\,\,\,\,(\forall j>1)=0$

Therefore,

$\lambda_{1}^{l}(\mathcal{M})=\max((\mathbf{S}_{l,m})_{\forall
m\in\theta_{1}(\mathcal{M})})\,\,\,\forall l\in\mathcal{L}$ (32)
$\lambda_{1}^{l}(\mathcal{M})=\pi_{1}^{l}(\mathcal{M})$

For all $l\in\mathcal{L}$ and $t=2$ we have,

$\lambda_{2}^{l}(\mathcal{M})=\max(\max(\mathbf{S}_{l,m})_{\forall
m\in\theta_{2}(\mathcal{M})},\max((\mathbf{S}_{l,m})_{\forall
m\in\theta_{1}(\mathcal{M})})*\rho,0)$ (33)

$\lambda_{2}^{l}(\mathcal{M})=\max(\max(\mathbf{S}_{l,m})_{\forall
m\in\theta_{2}(\mathcal{M})},\pi_{1}^{l}(\mathcal{M})*\rho)$ (34)
$\lambda_{2}^{l}(\mathcal{M})=\max((\mathbf{S}_{l,m})_{\forall
m\in\theta_{2}(\mathcal{M})},\rho\pi_{1}^{l}(\mathcal{M}))=\pi_{2}^{l}(\mathcal{M})$
(35)

For $l\in\mathcal{L}$ and $t=n$ suppose it is true that
$\lambda_{n}^{l}(\mathcal{M}))=\pi_{n}^{l}(\mathcal{M})$ and,

$\displaystyle\lambda_{n}^{l}(\mathcal{M})=\max($
$\displaystyle\max(\mathbf{S}_{l,m})_{\forall m\in\theta_{n}(\mathcal{M})},$
$\displaystyle\max((\mathbf{S}_{l,m})_{\forall
m\in\theta_{n-1}(\mathcal{M})})*\rho,\dots\max((\mathbf{S}_{l,m})_{\forall
m\in\theta_{1}(\mathcal{M})})*\rho^{n-1})$ (36)

For $l\in\mathcal{L}$ and $t=n+1$ we have,

$\displaystyle\lambda_{n+1}^{l}(\mathcal{M})=\max$
$\displaystyle(\max(\mathbf{S}_{l,m})_{\forall
m\in\theta_{n+1}(\mathcal{M})},$
$\displaystyle\max((\mathbf{S}_{l,m})_{\forall
m\in\theta_{n}(\mathcal{M})})*\rho,\dots\max((\mathbf{S}_{l,m})_{\forall
m\in\theta_{1}(\mathcal{M})})*\rho^{n})$ (37)

From Eq. 36,

$\lambda_{n+1}^{l}(\mathcal{M})=\max(\max(\mathbf{S}_{l,m})_{\forall
m\in\theta_{n+1}(\mathcal{M})},\rho*\lambda_{n}^{l}(\mathcal{M}))$ (38)

Hence,

$\lambda_{n+1}^{l}(\mathcal{M})=\max((\mathbf{S}_{l,m})_{\forall
m\in\theta_{n+1}(\mathcal{M})},\rho\pi_{n}^{l}(\mathcal{M}))=\pi_{n+1}^{l}(\mathcal{M})$
(39)

###  7.3 Prove of Theorem 2

Let $\theta_{t}$ be the set of locations that corresponds to sampled locations
in the sampling matrix $\bm{\Theta}$ for time stamp $t$. For $\pi_{0}=0$,
$\rho\in[0,1]$, the regressive facility location gain is defined by

$RFL(\mathcal{M})=\frac{\sum\limits_{t\in\mathcal{T}}\sum\limits_{l\in\mathcal{L}}\pi_{t}^{l}(\mathcal{M})}{L\times
T}$ (40)

where

$\pi_{t}^{l}(\mathcal{M})=\max((\mathbf{S}_{l,m})_{\forall
m\in\theta_{t}(\mathcal{M})},\rho\pi_{t-1}^{l}(\mathcal{M}))\,\,\,\forall
l\in\mathcal{L}$ (41)

####  7.3.1 Monotone non decreasing

A function f is submodular if
$\forall\,\,\,\mathcal{C}\subseteq\mathcal{D}\subseteq\mathcal{B}$

$f(\mathcal{C})\leq f(\mathcal{D})$

Proof:
Let $\mathcal{C}$ and $\mathcal{D}$ be the subset of buses following
$\forall\,\,\,\mathcal{C}\subseteq\mathcal{D}\subseteq\mathcal{B}$. We define
$\bm{\Theta}(\mathcal{C})$ and $\bm{\Theta}(\mathcal{D})$ as the sampling
matrix for the subset of buses $\mathcal{C}$ and $\mathcal{D}$.
Let $\theta_{t}$ be the set of locations that corresponds to 1 in the sampling
matrix $\bm{\Theta}(\mathcal{C})$ for the time stamp $t$ and Let
$\theta_{t}^{{}^{\prime}}$ be the set of locations that corresponds to 1 in
the sampling matrix $\bm{\Theta}(\mathcal{D})$ for the time stamp $t$.

$\theta_{t}=\\{i\\}\,\,\,\forall\,\,\,\bm{\Theta}_{i,t}(\mathcal{C})|_{=1}$
$\theta_{t}^{{}^{\prime}}=\\{i\\}\,\,\,\forall\,\,\,\bm{\Theta}_{i,t}(\mathcal{D})|_{=1}$

For all $l\in\mathcal{L}$ and $t=1$ we have

$\max\limits_{\forall
m\in\theta_{1}}(\mathbf{S}_{l,m})\leq\max\limits_{\forall
m\in\theta_{1}^{{}^{\prime}}}(\mathbf{S}_{l,m})$ (42)

For $t=2$ we have from using Eq. (22),

$\max((\mathbf{S}_{l,m})_{\forall m\in\theta_{2}},\rho\max\limits_{\forall
m\in\theta_{1}}(\mathbf{S}_{l,m}))\leq\max((\mathbf{S}_{l,m})_{\forall
m\in\theta_{2}^{{}^{\prime}}},\rho\max\limits_{\forall
m\in\theta_{1}^{{}^{\prime}}}(\mathbf{S}_{l,m}))$ (43)

for time stamp $t$ we have, from using Eq. (22)

$\displaystyle\max((\mathbf{S}_{l,m})_{\forall
m\in\theta_{t}},\rho\max\limits_{\forall
m\in\theta_{t-1}}(\mathbf{S}_{l,m}),\dots,\rho^{t-1}\max\limits_{\forall
m\in\theta_{1}}(\mathbf{S}_{l,m}))$
$\displaystyle\leq\max((\mathbf{S}_{l,m})_{\forall
m\in\theta_{t}^{{}^{\prime}}},\rho\max\limits_{\forall
m\in\theta_{t-1}^{{}^{\prime}}}(\mathbf{S}_{l,m}),\dots,\rho^{t-1}\max\limits_{\forall
m\in\theta_{1}^{{}^{\prime}}}(\mathbf{S}_{l,m}))$ (44)

Hence $f(\mathcal{C})\leq f(\mathcal{D})$

####  7.3.2 Submodular

A function f is submodular if
$\forall\,\,\,\mathcal{C}\subseteq\mathcal{D}\subseteq\mathcal{B}$ and
$b\in\mathcal{B}\backslash\mathcal{D}$

$f(\mathcal{C}\cup b)-f(\mathcal{C})\geq f(\mathcal{D}\cup b)-f(\mathcal{D})$

Consider a bus $b\in\mathcal{B}\backslash\mathcal{D}$, We add the bus $b$ in
both $\mathcal{C}$ and $\mathcal{D}$ and analyse the gain. Let the bus
traverse the set of location $G_{t}$ for a time stamp $t$.

$G_{t}=\\{i\\}\,\,\,\forall\,\,\,\bm{\Theta}_{i,t}(b)|_{=1}$

For all $l\in\mathcal{L}$ and $t=1\,\,\,\text{to}\,\,\,\mathcal{T}$ we have

$\displaystyle\pi_{t}^{l}(\mathcal{C}\cup b)-\pi_{t}^{l}(\mathcal{C})=$
$\displaystyle\max((\mathbf{S}_{l,m})_{\forall m\in\theta_{t}\cup
G_{t}},\rho\max\limits_{\forall m\in\theta_{t-1}\cup
G_{t-1}}(\mathbf{S}_{l,m}),\dots,\rho^{t-1}\max\limits_{\forall
m\in\theta_{1}\cup G_{1}}(\mathbf{S}_{l,m}))-$
$\displaystyle\max((\mathbf{S}_{l,m})_{\forall
m\in\theta_{t}},\rho\max\limits_{\forall
m\in\theta_{t-1}}(\mathbf{S}_{l,m}),\dots,\rho^{t-1}\max\limits_{\forall
m\in\theta_{1}}(\mathbf{S}_{l,m}))$ (45)

Let

$K=\max((\mathbf{S}_{l,m})_{\forall m\in\theta_{t}},\rho\max\limits_{\forall
m\in\theta_{t-1}}(\mathbf{S}_{l,m}),\dots,\rho^{t-1}\max\limits_{\forall
m\in\theta_{1}}(\mathbf{S}_{l,m}))$ $\displaystyle\pi_{t}^{l}(\mathcal{C}\cup
b)-\pi_{t}^{l}(\mathcal{C})=$ $\displaystyle\max((\mathbf{S}_{l,m})_{\forall
m\in\theta_{t}\cup G_{t}},\rho\max\limits_{\forall m\in\theta_{t-1}\cup
G_{t-1}}(\mathbf{S}_{l,m}),\dots,\rho^{t-1}\max\limits_{\forall
m\in\theta_{1}\cup G_{1}}(\mathbf{S}_{l,m}))-K$ (46)

$\displaystyle\pi_{t}^{l}(\mathcal{C}\cup b)-\pi_{t}^{l}(\mathcal{C})=$
$\displaystyle\max(0,(\mathbf{S}_{l,m})_{\forall m\in
G_{t}}-K,\rho\max\limits_{\forall m\in
G_{t-1}}(\mathbf{S}_{l,m})-K,\dots,\rho^{t-1}\max\limits_{\forall m\in
G_{1}}(\mathbf{S}_{l,m})-K)$ (47)

let

$K^{\prime}=\max((\mathbf{S}_{l,m})_{\forall
m\in\theta_{t}^{{}^{\prime}}},\rho\max\limits_{\forall
m\in\theta_{t-1}^{{}^{\prime}}}(\mathbf{S}_{l,m}),\dots,\rho^{t-1}\max\limits_{\forall
m\in\theta_{1}^{{}^{\prime}}}(\mathbf{S}_{l,m}))$

From Eq. 44 since $K\leq K^{\prime}$ Therefore,

$\displaystyle\pi_{t}^{l}(\mathcal{C}\cup b)-\pi_{t}^{l}(\mathcal{C})\geq$
$\displaystyle\max(0,(\mathbf{S}_{l,m})_{\forall m\in
G_{t}}-K^{\prime},\rho\max\limits_{\forall m\in
G_{t-1}}(\mathbf{S}_{l,m})-K^{\prime},\dots,\rho^{t-1}\max\limits_{\forall
m\in G_{1}}(\mathbf{S}_{l,m})-K^{\prime})$ (48)
$\displaystyle\pi_{t}^{l}(\mathcal{C}\cup b)-\pi_{t}^{l}(\mathcal{C})\geq$
$\displaystyle\max((\mathbf{S}_{l,m})_{\forall m\in\theta_{t}\cup
G_{t}},\rho\max\limits_{\forall m\in\theta_{t-1}\cup
G_{t-1}}(\mathbf{S}_{l,m}),\dots,\rho^{t-1}\max\limits_{\forall
m\in\theta_{1}\cup G_{1}}(\mathbf{S}_{l,m}))-K^{\prime}$ (49)

$\displaystyle\pi_{t}^{l}(\mathcal{C}\cup b)-\pi_{t}^{l}(\mathcal{C})\geq$
$\displaystyle\max((\mathbf{S}_{l,m})_{\forall m\in\theta_{t}^{\prime}\cup
G_{t}},\rho\max\limits_{\forall m\in\theta_{t-1}^{\prime}\cup
G_{t-1}}(\mathbf{S}_{l,m}),\dots,\rho^{t-1}\max\limits_{\forall
m\in\theta_{1}^{\prime}\cup G_{1}}(\mathbf{S}_{l,m}))$
$\displaystyle-\max((\mathbf{S}_{l,m})_{\forall
m\in\theta_{t}^{\prime}},\rho\max\limits_{\forall
m\in\theta_{t-1}^{\prime}}(\mathbf{S}_{l,m}),\dots,\rho^{t-1}\max\limits_{\forall
m\in\theta_{1}^{\prime}}(\mathbf{S}_{l,m}))$ (50)
$\displaystyle\pi_{t}^{l}(\mathcal{C}\cup
b)-\pi_{t}^{l}(\mathcal{C})\geq\pi_{t}^{l}(\mathcal{D}\cup
b)-\pi_{t}^{l}(\mathcal{D})$ (51)

From Eqs (45-51) for all $l\in\mathcal{L}$ and
$t=1\,\,\,\text{to}\,\,\,\mathcal{T}$ we have,

$f(\mathcal{C}\cup b)-f(\mathcal{C})\geq f(\mathcal{D}\cup b)-f(\mathcal{D})$

###  7.4 Creating Dense Maps from Sampled Data

We use Variational Bayesian Matrix Completion (Cold start) and Variational
Bayesian Subspace Filtering (Cold start) as the extended matrix completion
frameworks that handles the cold start location data imputation to evaluate
the performance of the dense AQ maps. We impute the missing data using
VBMC(CS) where we optimize the Eqs. (52, 54). We also impute the missing data
using VBSF(CS) where we optimize the Eqs. (52, 53, 54). The matrix completion
framework (VBMC) proposed in paper [20] imposes a low rank structure on the
data to impute the missing data as:

$\displaystyle\mathcal{L}_{1}=\min_{\mathbf{A},\mathbf{B}}||\mathbf{P}_{\Omega}(\mathbf{Y}-\mathbf{A}\mathbf{B}^{T})||_{{F}}$
(52)

where $\mathbf{A}\in\mathcal{R}^{L\times r}$ and
$\mathbf{B}\in\mathcal{R}^{T\times r}$ and $r=$rank$(\mathbf{Y})<<min(L,T)$
Further, VBSF [5] add a regularization on the matrix $\mathbf{B}$ to
incorporate the temporal evolution in addition to the low rankness (Eq. 52) as

$\displaystyle\mathcal{R}(\mathbf{B})=\sum_{i=1}^{T}||\mathbf{b}_{i}-\mathbf{F}\mathbf{b}_{i-1}||$
(53)

However, these framework does not incorporate for the cold start locations.
Incorporating the similarity matrix $\mathbf{G}$ along with the low rank
matrix completion framework can tackle the cold start problem [23, 24].

$\displaystyle\mathcal{L}_{2}=\min_{\mathbf{A},\mathbf{C}}||(\mathbf{G}-\mathbf{A}\mathbf{C}^{T})||_{{F}}$
(54)

Refer the update equations in the papers [5, 20]. We showed the changed
updated equations below:

The update for $i^{th}$ column of $\mathbf{A}$ for the cold start matrix
completion is as follows:

$\displaystyle\bm{\Xi}^{\mathbf{A}}_{i}$
$\displaystyle=\left(\hat{\gamma}_{i}\mathbf{I}_{r}+\hat{\beta}\sum_{\tau\in\Omega^{\prime}_{i}}(\bm{\mu}_{\tau}^{\mathbf{B}}(\bm{\mu}_{\tau}^{\mathbf{B}})^{T}+\bm{\Xi}_{\tau,\tau}^{\mathbf{B}})+\hat{\beta_{1}}(\bm{\mu}^{\mathbf{C}}(\bm{\mu}^{\mathbf{C}})^{T}+\bm{\Xi}^{\mathbf{C}})\right)^{-1}$
(55) $\displaystyle\bm{\mu}^{\mathbf{A}}_{i}$
$\displaystyle=\bm{\Xi}^{\mathbf{A}}_{i}(\hat{\beta}\sum_{\tau\in\Omega^{\prime}_{i}}\bm{\mu}^{\mathbf{B}}_{\tau}y_{i\tau}+\hat{\beta_{1}}\bm{\mu}^{\mathbf{C}}g_{i})$
(56)

The update for $\beta_{1}$ is as follows:

$\hat{\beta_{1}}=\frac{p\,L^{2}}{\parallel\mathbf{G}-\mathbf{A}\mathbf{C}^{T}\parallel^{2}_{F}}$
(57)

The updates Eqs. for $\mathbf{C}$ is as follows:

$\displaystyle\Gamma$ $\displaystyle=diag(\gamma),$
$\displaystyle\bm{\Xi}^{\bf C}\tt$ $\displaystyle=(<\beta_{1}><\bf
A^{T}A\tt>+\Gamma)^{-1}$ (58) $\displaystyle<\mathbf{C}>$
$\displaystyle=<\beta_{1}>\mathbf{G}\it<\mathbf{A}>\bm{\Xi}^{\bf C},$ (59)