Decentralized Data Fusion and Active Sensing with Mobile Sensors for Modeling and Predicting Spatiotemporal Traffic Phenomena
\hspace{-12mm}{Jie Chen$^{\dag}$
Kian Hsiang Low$^{\dag}$
Colin Keng-Yan Tan$^{\dag}$,}
\hspace{-12mm}{Ali Oran$^{\S}$}
\hspace{-12mm}Dept. Computer Science
\hspace{-12mm}Singapore-MIT Alliance for Research and Technology$^{\S}$
\hspace{-12mm}Republic of Singapore \And \hspace{12mm} {Patrick Jaillet$^{\ddag}$
John Dolan$^{\P}$
Gaurav Sukhatme$^{\ast}$}
\hspace{12mm}Dept. Electrical Engineering and Computer Science
MIT$^{\ddag}$
\hspace{12mm}Dept. Computer Science
\hspace{12mm}USA %\And %{Coauthor}
%Affiliation
%Address
%(if needed)
\vspace{-5mm} The problem of modeling and predicting spatiotemporal traffic phenomena over an urban road network is important to many traffic applications such as detecting and forecasting congestion hotspots. This paper presents a decentralized data fusion and active sensing (D$^2$FAS) algorithm for mobile sensors to actively explore the road network to gather and assimilate the most informative data for predicting the traffic phenomenon. We analyze the time and communication complexity of D$^2$FAS and demonstrate that it can scale well with a large number of observations and sensors. We provide a theoretical guarantee on its predictive performance to be equivalent to that of a sophisticated centralized sparse approximation for the Gaussian process (GP) model: The computation of such a sparse approximate GP model can thus be parallelized and distributed among the mobile sensors (in a Google-like MapReduce paradigm), thereby achieving efficient and scalable prediction. We also theoretically guarantee its active sensing performance that improves under various practical environmental conditions. Empirical evaluation on real-world urban road network data shows that our D$^2$FAS algorithm is significantly more time-efficient and scalable than state-of-the-art centralized algorithms while achieving comparable predictive performance.\vspace{-3mm}
uai2012.tex
% Sample LaTeX file for creating a paper in the Morgan Kaufmannn two 
% column, 8 1/2 by 11 inch proceedings format. 

\documentclass[]{article}
\usepackage{proceed2e}
\usepackage{times}
\usepackage{natbib}
%\usepackage{named}
\usepackage[ruled,vlined]{algorithm2e}
\usepackage{amsmath}
\usepackage{amssymb}
%\usepackage{dsfont}
\usepackage{graphicx}
\usepackage{wrapfig}

\title{Decentralized Data Fusion and Active Sensing with Mobile Sensors for Modeling and Predicting Spatiotemporal Traffic Phenomena} 

\author{\hspace{-12mm}{\bf Jie Chen$^{\dag}$, Kian Hsiang Low$^{\dag}$, Colin Keng-Yan Tan$^{\dag}$,}\\
\hspace{-12mm}{\bf Ali Oran$^{\S}$}\\
       \hspace{-12mm}Dept. Computer Science, National University of Singapore$^\dag$\\
       \hspace{-12mm}Singapore-MIT Alliance for Research and Technology$^{\S}$\\
       \hspace{-12mm}Republic of Singapore
\And \hspace{12mm}
{\bf Patrick Jaillet$^{\ddag}$, John Dolan$^{\P}$, Gaurav Sukhatme$^{\ast}$}\\  
       \hspace{12mm}Dept. Electrical Engineering and Computer Science, MIT$^{\ddag}$\\
       \hspace{12mm}The Robotics Institute, Carnegie Mellon University$^{\P}$\\
       \hspace{12mm}Dept. Computer Science, University of Southern California$^{\ast}$\\
       \hspace{12mm}USA            
%\And 
%{\bf Coauthor}   \\ 
%Affiliation \\          
%Address    \\           
%(if needed)\\ 
} 

\newcounter{sol} 
\setcounter{sol}{1} % with solutions
%\setcounter{sol}{0} % without solutions
 
\begin{document} 
\newcommand{\set}[1]{{\mathcal #1}}
\newcommand{\squishlisttwo}{
 \begin{list}{$\bullet$}
  { \setlength{\itemsep}{0pt}
     \setlength{\parsep}{0pt}
    \setlength{\topsep}{0pt}
    \setlength{\partopsep}{0pt}
    \setlength{\leftmargin}{1em}
    \setlength{\labelwidth}{1.5em}
    \setlength{\labelsep}{0.5em} } }

\newcommand{\squishend}{
  \end{list}  }

\newtheorem{defi}{Definition}
\newtheorem{thm}{Theorem}
 
\maketitle 
 
\begin{abstract}
\vspace{-5mm}
The problem of modeling and predicting spatiotemporal traffic phenomena over an urban road network is important to many traffic applications such as detecting and forecasting congestion hotspots. This paper presents a decentralized data fusion and active sensing (D$^2$FAS) algorithm for mobile sensors to actively explore the road network to gather and assimilate the most informative data for predicting the traffic phenomenon. 
We analyze the time and communication complexity of D$^2$FAS and demonstrate that it can scale well with a large number of observations and sensors.
We provide a theoretical guarantee on its predictive performance to be equivalent to that of a sophisticated centralized sparse approximation for the Gaussian process (GP) model: The computation of such a sparse approximate GP model can thus be parallelized and distributed among the mobile sensors (in a Google-like MapReduce paradigm), thereby achieving efficient and scalable prediction.
We also theoretically guarantee its active sensing performance that improves under various practical environmental conditions.
Empirical evaluation on real-world urban road network data shows that our D$^2$FAS algorithm
is significantly more time-efficient and scalable than state-of-the-art centralized algorithms while achieving comparable predictive performance.\vspace{-3mm}
\end{abstract}

\section{Introduction}
\label{sec:intro}
\vspace{-3mm}
%
Knowing and understanding the traffic conditions and phenomena over road networks has become increasingly important to the goal of achieving smooth-flowing, congestion-free traffic, especially in densely-populated urban cities.
According to a $2011$ urban mobility report \citep{Schrank11}, traffic congestion in the USA has caused $1.9$ billion gallons of extra fuel, $4.8$ billion hours of travel delay, and $\$101$ billion of delay and fuel cost. 
Such huge resource wastage can potentially be mitigated if the spatiotemporally varying traffic phenomena (e.g., speeds and travel times along road segments) 
%over an urban road network 
are predicted accurately enough in real time to detect and forecast the congestion hotspots; network-level (e.g., ramp metering, road pricing) and user-level (e.g., route replanning) measures can then be taken to relieve this congestion, so as to improve the overall efficiency of road networks.

In practice, it is non-trivial to achieve real-time, accurate prediction of a spatiotemporally varying traffic phenomenon because the quantity of sensors that can be deployed to observe an entire road network is cost-constrained.
Traditionally, static sensors such as loop detectors \citep{krause08,wang05} are placed at designated locations in a road network to collect data for %learning a spatiotemporal predictive model of 
predicting the traffic phenomenon.
% The massive amount of wasted resources indicates that substantial gain can be potentially obtained if the accurate and reliable real-time traffic phenomena (e.g. velocity) is available. With the knowledge of traffic phenomena, on the system level, congestion can be controlled by ramp metering or dynamic speed limits; on the individual level, travel path or travel time can be planned avoiding congestion. 
% Consequently, the overall efficiency of the road network is improved. 
% Historically, the traffic data for phenomena modeling and prediction rely on dedicate sensing infrastructure (loop detectors, radars and cameras etc.).
However, they provide sparse coverage (i.e., many road segments are not observed, thus leading to data sparsity), incur high installation and maintenance costs, and cannot reposition by themselves in response to changes in the traffic phenomenon.
%static positions disabling optimized or adaptive redeployment when traffic phenomena changed.
Low-cost GPS technology allows the collection of traffic data using passive mobile probes \citep{Bayen10a} (e.g., taxis/cabs). 
%krause08,
Unlike static sensors, they can directly measure the travel times along road segments. But, they provide fairly sparse coverage due to low GPS sampling frequency (i.e., often imposed by taxi/cab companies) and no control over their routes, incur high initial implementation cost, pose privacy issues, and produce highly-varying speeds and travel times while traversing the same road segment due to inconsistent driving behaviors. A critical mass of probes is needed on each road segment to ease the severity of the last drawback \citep{Srinivasan96} but is often hard to achieve on non-highway segments due to sparse coverage.
% high initial implementation costs.
In contrast, we propose the use of active mobile probes \citep{turner98} to overcome the limitations of static and passive mobile probes. In particular, they can be directed to explore any segments of a road network to gather traffic data at a desired GPS sampling rate while enforcing consistent driving behavior.
%a small number of 
%the active mobile probes 
% the most informative data for %learning a spatiotemporal predictive model of 
%modeling and predicting the traffic phenomenon.
% Recently, many works are devoted to mobile traffic sensing methods which can overcome the limitations of dedicated sensing infrastructure. However, these methods rely on mobile sensors which bring in new challenges.  The mobility of sensors enables more road segments (highway and arterial) be covered, which indicates a more dense road network with complex topology causing failure of existing models.  The mobile sensors need to jointly determine the most 'informativ\mathcal{E} paths to gather the traffic data.  The gathered traffic data are distributed among mobile sensors, therefore, have to be assimilated so that can be transformed into traffic phenomena.  Since the number of mobile sensors possibly scales up to a large amount, scalability and real time property are worthy of crucial attention. 

How then do the mobile probes/sensors actively explore a road network to gather and assimilate the most informative observations for predicting the traffic phenomenon? There are three key issues surrounding this problem, which will be discussed together with the related works:

{\bf Models for predicting spatiotemporal traffic phenomena.}
%{\bf Spatiotemporal traffic prediction models.}
%spatia-temporal models our approach GP-GK
%In existing literatures, spatiotemporal models are frequently used to model traffic phenomena, since traditional temporal model can only capture the evolution of the traffic phenomena but cannot deal with missing data on segment of road network. Therefore, spatial correlation is taken into account in certain manner. 
The spatiotemporal correlation structure of a traffic phenomenon can be exploited to predict the traffic conditions of any unobserved road segment at any time using the observations taken along the sensors' paths. To achieve this, existing Bayesian filtering frameworks \citep{chen11,wang05,Bayen10a}
%mihaylova07a,
utilize various handcrafted parametric models predicting traffic flow along a highway stretch
%utilize various handcrafted parametric traffic flow models of a highway stretch
that only correlate adjacent segments of the highway.
So, their predictive performance will be compromised when the current observations are sparse and/or the actual spatial correlation spans multiple segments.
%they may not fit the true underlying correlation structure well, thus compromising their predictive performance, especially with limited observations.
Their strong Markov assumption further exacerbates this problem.
It is also not shown how these models can be generalized to work for arbitrary road network topologies and more complex correlation structure.
%assume observations to be conditionally independent across highway segments
%traffic flow model
%Another possible model is the macroscopic traffic flow model combining with recursive Bayesian filter as prediction model, in which particle filter incurs intensive time and storage complexity \citep{chen11,mihaylova07a} and extended Kalman filters \citep{wang05} cannot track multi-model distributions.
Existing multivariate parametric traffic prediction models \citep{kamarianakis03,min11} do not quantify uncertainty estimates of the predictions and impose rigid spatial locality assumptions that do not adapt to the true underlying correlation structure.
%require storing and estimating a large number of parameters (i.e., quadratic in the number of road segments), 
%parametric
%integrate the spatial correlation as the component-wise product of parameter matrix and weight matrix reflecting configuration of road network. However, these models do not exploit the features (e.g. length, width, speed limit, direction etc.) of segments for prediction.  These models tend to be computational intensive as the number of spatial correlation parameters is quadratic in number of road segments.
%In \citep{blandin09b}, additional features are represented as a set of covariants, and the most relevant features can be selected via optimization. However, this method requires the entry time known on each segment and the spatial correlation is not utilized.  A limitation of the above methods is that they cannot provide uncertainty measure.  
%bayes network
%\citep{sun06} incorporates spatial correlation by modeling the road network as a Bayesian network. The Markov assumption simplifies the prediction model by reducing the number of cause nodes. However, the sparse traffic data do not favor a high dimension predictive distribution.  

In contrast, we assume the traffic phenomenon over an urban road network (i.e., comprising the full range of road types like highways, arterials, slip roads) to be realized from a rich class of Bayesian non-parametric models called the \emph{Gaussian process} (GP) (Section~\ref{sec:gpgk}) that can formally characterize its spatiotemporal correlation structure and be refined with growing number of observations.
% \citep{Teh10}. 
More importantly, GP can provide formal measures of predictive uncertainty (e.g., based on variance or entropy criterion) for directing the sensors to explore highly uncertain areas of the road network. 
\cite{krause08} used GP to represent the traffic phenomenon over a network of only highways and defined the correlation of speeds between highway segments to depend only on the geodesic (i.e., shortest path) distance of these segments with respect to the network topology; their features are not considered.
%like the work of \citep{krause08}, 
\cite{Kersting09} maintained a mixture of two independent GPs for flow prediction such that the correlation structure of one GP utilized road segment features while that of the other GP depended on manually specified relations (instead of geodesic distance) between segments with respect to an undirected network topology.
Different from the above works, we propose a relational GP
% \citep{Chu07} 
whose correlation structure exploits the geodesic distance between segments based on the topology of a directed road network with vertices denoting road segments and edges indicating adjacent segments weighted by dissimilarity of their features, hence tightly integrating the features and relational information.
%embeds road segment features information as edge costs in a weighted directed graph representing the road network
% is based on a tight coupling of both the features and topology information by exploiting a weighted directed graph with feature-dependent edge costs and geodesic distance defined relative to its topology.
%Different from work of \cite{krause08}, we further improve the correlation structure of GP by enabling it to exploit road segment features 
%(e.g., length, number of lanes, direction, speed limit) 
%for differentiating road types, which is not found in the works described above.
%It is desirable to develop prediction models that can exploit , topology of road network, and the road features.
% uncertainty measures for confidence estimates and to direct robots
%\citep{krause08} uses Gaussian process (GP) to model the spatiotemporal phenomena, in which the correlation of velocity on road segments depends on the geodesic (i.e. shortest path) distance of these segments with respect to the road network topology. This model is close to the one we use, while we augment it by taking into account the features of road segments.  Based on the GP model, our approach relates to the existing techniques in two aspects:
%kernel framework
%closest model:GP-GK

{\bf Data fusion.} 
The observations are gathered distributedly by each sensor along its path in the road network and have to be assimilated in order to predict the traffic phenomenon. Since a large number of observations are expected to be collected, a centralized approach to GP prediction cannot be performed in real time due to its cubic time complexity.

To resolve this, we propose a decentralized data fusion approach to efficient and scalable approximate GP prediction (Section~\ref{sec:ddf}).
%computational load for each sensor decreases with increasing number of robots.
% that, in particular, exploits the spatiotemporal correlation structure of a traffic phenomenon. 
Existing decentralized and distributed Bayesian filtering frameworks for addressing non-traffic related problems \citep{Chung04,Coates04,saber05a,Thrun03,Durrant-Whyte03b}
will face the same difficulties as their centralized counterparts described above if applied to predicting traffic phenomena, thus resulting in loss of predictive performance.
% adopt parametric models and strong Markov assumption that may not approximate the true underlying correlation structure well.
%mihaylova12,
%data loss during communication.
Distributed regression algorithms \citep{guestrin04,Guestrin04b} for static sensor networks gain efficiency from spatial locality assumptions, which cannot be exploited by mobile sensors whose paths are not constrained by locality.
%are easily violated since each mobile sensor's path does not have to be constrained to a locality.
\cite{cortes09} proposed a distributed data fusion approach to approximate GP prediction based on an iterative Jacobi overrelaxation algorithm, which incurs some critical limitations: (a) the past observations taken along the sensors' paths are assumed to be uncorrelated, which greatly undermines its predictive performance when they are in fact correlated and/or the current observations are sparse;
%thus greatly undermining the predictive performance; 
(b) when the number of sensors grows large, it converges very slowly; 
(c) it assumes that the range of positive correlation has to be bounded by some factor of the communication range. Our proposed decentralized algorithm does not suffer from these limitations and can be computed exactly with efficient time bounds.
%are violated since the mobile sensors' paths are not localized.
% With multiple sensors, decentralized data fusion algorithms are to assimilate the distributed local data into a consistent view so that the mobile sensors can make joint decision consistently.  In distributed particle filters, the consistency can be achieved by either transmitting particle values and their weights between the sensors \citep{sheng05} or communicating a parametric approximation of the conditional density \citep{mihaylova12}. Our approach is more close to the later one. Similar approach is also adopted in distributed Kalman filter \citep{saber05a}, in which the data are fused in term of different natures: weighted measurement and inverse-covariance matrices.  In our approach, local data are summarized into similar structure of summarized measurements and covariance-like matrices, and exchanged to form a consistent view.   Another close approach is the kernel linear regression \citep{guestrin04} where the local correlation in data expressed by coefficients given a set of basis functions. In our approach, the support set assisting local data summarization consists of a small set of segments which is similar to the basis functions in kernel linear regression. 

{\bf Active sensing.} 
The sensors have to coordinate to actively gather the most informative observations for minimizing the uncertainty of modeling and predicting the traffic phenomenon. Existing centralized \citep{LowAAMAS08,LowICAPS09,LowAAMAS11} and decentralized \citep{,LowAAMAS12,stranders09} active sensing algorithms scale poorly with a large number of observations and sensors.
%krause08,Guestrin08,
We propose a partially decentralized active sensing algorithm that overcomes these issues of scalability (Section~\ref{sec:dadgf}).
%\vspace{1mm} 
%\citep{stranders09} presents a decentralized coordination algorithm for mobile sensors performing active sensing based on GP model.  However, this approach never improves the efficiency of computation w.r.t size of observed data, which is prohibitively expensive due to the large number of joint entropy evaluation.  Moreover, the joint action pruning algorithm, in the worst case, is still exponential in the number of agents and length of planning horizon. In contrast, our approach relies on the support set summarizing the local data so that the joint entropy evaluation does not depend on the whole data directly. Moreover, the joint entropy computation is fully decentralized due to conditional independence property.
%test on simulated data, not provide computation and communication
%complexity
%\citep{krause08} aims at predicting traffic phenomena with static sensors and passive mobile sensors. Unlike our active sensing approach, traffic data for phenomena modeling are gathered from the passive sensors whose importance is quantified by variance reduction.  
%\citep{turner98} provides the in-depth comparison between approaches using passive and active mobile sensors. Our approach can also be flexibly applied to passive mobile sensors. \\

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%contribution
This paper presents a novel \emph{Decentralized Data Fusion and Active Sensing} (D$^2$FAS) algorithm (Sections~\ref{sec:ddf} and~\ref{sec:dadgf}) for sampling spatiotemporally varying environmental phenomena with mobile sensors. 
Note that the decentralized data fusion component of D$^2$FAS can also be used for static and passive mobile sensors.
The practical applicability of D$^2$FAS is not restricted to traffic monitoring; it can be used in other environmental sensing applications such as mineral prospecting \citep{LowICRA07}, monitoring of ocean and freshwater phenomena \citep{LowSPIE09,LowAeroconf10,LowIPSN09,LowAAMAS11,LowAAMAS12} (e.g., plankton bloom, anoxic zones), forest ecosystems, pollution (e.g., oil spill), or contamination (e.g., radiation leak).
%, which can be applied to mobile sensors to actively gather traffic data from road network and assimilate distributed data into a global view so that joint exploring strategy can be determined consistently. The traffic phenomena is modeled by Gaussian process (GP) with graph-based kernel (Section \ref{sec:gpgk}), a spatiotemporal model which is suitable for traffic modeling and prediction. Additionally, approximation of GP is introduced to balance time efficiency and prediction accuracy.  
The specific contributions of this paper include:\vspace{-0.6mm}
%\begin{itemize}
\squishlisttwo
\item Analyzing the time and communication overheads of D$^2$FAS (Section~\ref{overhead}): We prove that D$^2$FAS can scale better than existing state-of-the-art centralized algorithms with a large number of observations and sensors;
%outperforms the existing state-of-art centralized approaches in terms of efficiency and scalability.
\item Theoretically guaranteeing the predictive performance of the decentralized data fusion component of D$^2$FAS to be equivalent to that of 
%a sophisticated centralized approximate GP prediction model 
a sophisticated centralized sparse approximation for the GP model (Section~\ref{sec:ddf}): The computation of such a sparse approximate GP model can thus be parallelized and distributed among the mobile sensors (in a Google-like MapReduce paradigm), thereby achieving efficient and scalable prediction;
\item Theoretically guaranteeing the performance of the partially decentralized active sensing component of D$^2$FAS, from which various practical environmental conditions can be established to improve its performance;
%\item Improving the correlation structure of GP model by enabling it to exploit road segment features (e.g., length, number of lanes, direction, and speed limit) and the road network topology (Section~\ref{sec:gk});
\item Developing a relational GP model whose correlation structure can exploit both the road segment features
% (e.g., length, number of lanes, direction, and speed limit) 
and road network topology information (Section~\ref{sec:gk});
\item Empirically evaluating the predictive performance, time efficiency, and scalability of the D$^2$FAS algorithm on a real-world traffic phenomenon (i.e., speeds of road segments) dataset over an urban road network (Section~\ref{sec:exp}): D$^2$FAS is more time-efficient and scales significantly better with increasing number of observations and sensors while achieving predictive performance close to that of existing state-of-the-art centralized algorithms. \vspace{-1mm}
\squishend
%\end{itemize}
%
\section{Relational Gaussian Process Regression}
\label{sec:gpgk}
\vspace{-3mm}
%Gaussian Process Regression
%Gaussian process (GP) \citep{rasmussen06} is an effective spatiotemporal model which can estimate traffic phenomena on the unobserved segments given a set of observed segments, meanwhile quantify the uncertainty of the prediction. 
%Additionally, a well-known limitation of GP is that it requires cubic time in the size of observations to do prediction.  In this section, we will briefly introduce GP model, the graph-based kernel model over road network and an approximation of GP based on subset of data.
% Different from normal GP, we need to take into account road segment features and network topology.
%
%\subsection{Gaussian Process}
%\label{sec:gp} 
%
The \emph{Gaussian process} (GP) can be used to model a spatiotemporal traffic phenomenon over a road network as follows: The traffic phenomenon is defined to vary as a realization of a GP.
Let $V$ be a set of road segments representing the domain of the road network such that each road segment $s\in V$ is specified by a $p$-dimensional vector of features and is associated with a realized (random) measurement $z_s$ ($Z_s$) of the traffic condition such as speed if $s$ is observed (unobserved). Let $\{Z_s\}_{s \in V}$ denote a GP, that is, every finite subset of $\{Z_s\}_{s \in V}$ follows a multivariate Gaussian distribution \citep{rasmussen06}. Then, the GP is fully specified by its \emph{prior} mean $\mu_{s} \triangleq \mathbb{E}[Z_s]$ and covariance $\sigma_{s s'} \triangleq \mbox{cov}[Z_s, Z_{s'}]$ for all $s, s' \in V$. 
In particular, we will describe in Section~\ref{sec:gk} how the covariance $\sigma_{s s'}$ for modeling the correlation of measurements between all pairs of segments $s, s' \in V$ can be designed to exploit the road segment features and the road network topology.
%However, in traffic domain, certain amount of correlation is hidden due to topology of road network, while the existing kernel functions cannot model it well.  
%GP models the domain as an infinite multivariate Gaussian distribution:
%
%\begin{equation}
%Z_{V}\sim \mathcal{N}( M,\Sigma)
%\label{eq:gp}
%\end{equation}
%
%where $Z_{V}$ are the random variables over domain $V$.  Here GP model is fully specified by a mean function $ M(.)$ and a kernel function (also called covariance function) $\Sigma(.,.)$ which are respectively defined as
%
%\begin{equation}
% M:V
%\mapsto
%\mathbb{R},
%\label{eq:mean}
%\end{equation}
%
%\begin{equation}
%\Sigma:V\times V\mapsto\mathbb{R}.
%\label{eq:cov}
%\end{equation}
%
%Since the measurements on the segments are noisy, we define $f_s=Z_s+\epsilon$ where $s\in V$ and $\epsilon\sim\mathcal{N}(0,\sigma_n^2)$.  

A chief capability of the GP model is that of performing probabilistic regression:
Given a set $D\subset V$ of observed road segments and a column vector $z_D$ of corresponding measurements, the joint distribution of the measurements 
%$Z_Y$ 
at any set $Y\subseteq V\setminus D$ of unobserved road segments remains Gaussian with the following \emph{posterior} mean vector and covariance matrix\vspace{-0mm}
%
\begin{equation}
\displaystyle\mu_{Y|D}\triangleq \mu_Y + \Sigma_{YD}\Sigma_{DD}^{-1}( z_D - \mu_D )
\label{pmu}\vspace{-0mm}
\end{equation}
\begin{equation}
\displaystyle\Sigma_{YY|D}\triangleq\Sigma_{YY}-\Sigma_{YD}\Sigma_{DD}^{-1}\Sigma_{DY}
\label{pcov}\vspace{-0mm}
\end{equation}
%
where ${\mu}_Y$ (${\mu}_D$) is a column vector with mean components $\mu_{s}$ for all $s\in Y$ ($s\in D$),
$\Sigma_{YD}$ ($\Sigma_{DD}$) is a covariance matrix with covariance components $\sigma_{s s'}$ for all $s\in Y, s^\prime\in D$ ($s, s^\prime\in D$), and $\Sigma_{DY}$ is the transpose of $\Sigma_{YD}$. The posterior mean vector $\mu_{Y|D}$ (\ref{pmu}) is used to predict the measurements at any set $Y$ of unobserved road segments.
The posterior covariance matrix $\Sigma_{YY|D}$ (\ref{pcov}), which is independent of the measurements $z_D$,
can be processed in two ways to quantify the uncertainty of these predictions: (a) the trace of $\Sigma_{YY|D}$ yields the sum of posterior variances $\Sigma_{ss|D}$ over all $s\in Y$; (b) the determinant of $\Sigma_{YY|D}$ is used in calculating the Gaussian posterior joint entropy
\vspace{-2.5mm}
%
\begin{equation}
  \mathbb{H}[Z_Y|Z_D]\triangleq
  \frac{1}{2}\log\hspace{-0.5mm}\left(2\pi e\right)^{\left|Y\right|}\left|\Sigma_{YY|D}\right| \ .
  \label{entropy}\vspace{-2.5mm}
\end{equation}
%
In contrast to the first measure of uncertainty that assumes conditional independence between measurements in the set $Y$ of unobserved road segments, the entropy-based measure (\ref{entropy}) accounts for their correlation, thereby not overestimating their uncertainty.
Hence, we will focus on using the entropy-based measure of uncertainty in this paper.\vspace{-2mm}
% we can perform GP regression on an unobserved set $Y\subset V\setminus D$.  
%conditional density
%An important property of GP is that any finite collection of random variables follows a Gaussian distribution. 
% Given two sets $A,B\subset V$, we denote $\Sigma_{A,B}$ as their kernel matrix whose elements are evaluated by (\ref{eq:cov}).  
%GP regression can be performed by predicting the measurements of unobserved road segments using the posterior mean (\ref{pmu}) and the variance of these predictions correspond to the diagonal components of the corresponding posterior covariance (\ref{pcov}).
%A larger entropy indicates that the prediction is of higher uncertainty, or in information theory, the set $Y$ is more informative. Mobile sensors can exploit this feature to measure unobserved segments with highest uncertainty so that the uncertainty reduction in the domain is the maximum.
%
\subsection{Graph-Based Kernel}
\label{sec:gk}\vspace{-3mm}
If the observations are noisy (i.e., by assuming additive independent identically distributed Gaussian noise with variance $\sigma^2_n$), then their prior covariance $\sigma_{s s'}$ can be expressed as
$
\sigma_{s s'} = k(s,s') + \sigma^2_n \delta_{s s'}
$
%for all $s, s' \in V$ 
where $\delta_{s s'}$ is a Kronecker delta that is $1$ if $s = s'$ and $0$ otherwise, and $k$ is a kernel function measuring the pairwise ``similarity'' of road segments.
%which will be formally defined next.
%In GP model, the most important component is kernel function (\ref{eq:cov}) which captures the correlation between segment pair.
%Traditionally, the kernel function $k$ measures the ``similarity'' between two feature vectors.
%However, 
For a traffic phenomenon (e.g., road speeds),
the correlation of measurements between pairs of road segments depends not only on their features (e.g., length, number of lanes, speed limit, direction) but also
%the kernel value is evaluated from kernel function with feature vectors of two segments as the inputs.  
%in terms of traffic phenomena (specifically velocity on road segments), the correlation of segment pair not only depends on the information provided
%in feature vector    but also
the road network topology.  
So, the kernel function is defined to exploit both the features and topology information, which will be described next.\vspace{-1.5mm}  
\begin{defi}[Road Network]
Let the road network be represented as a weighted directed graph $G\triangleq (V,E,m)$ that consists of\vspace{-1mm}
\squishlisttwo
\item a set $V$ of vertices denoting the domain of all possible road segments,
\item a set $E\subseteq V\times V$ of edges where there is an edge $(s,s')$ from $s\in V$ to $s'\in V$ iff the end of segment $s$ connects to the start of segment $s'$ in the road network, and
\item a weight function $m: E\rightarrow \mathbb{R}^{+}$ measuring 
the standardized Manhattan distance \citep{borg05}
%\vspace{-3mm}
%\begin{equation*}
$m((s,s'))\triangleq\sum_{i=1}^p \left|[s]_i-[s']_i\right|\hspace{-0.5mm}/r_i$
%$w((s,s'))=\sum_{i=1}^p |s_i-s'_i|/r_i$
%\label{eq:w}
%\end{equation*}
 of each edge $(s,s')$
where $[s]_i$ ($[s']_i$) is the $i$-th component of the feature vector specifying road segment $s$ ($s'$), and $r_i$ is the range of the $i$-th feature. 
The weight function $m$ serves as a dissimilarity measure between adjacent road segments.\vspace{-1.5mm}
\squishend
\end{defi} 
%the modified kernel function be denoted as $k(s,s';G)$ where 
%
%\begin{equation}
 % G=(V,E)
%\label{eq:g}
%\end{equation}
%
%where $V$ is the domain of all possible road segments, 
%indicates the connectivity between segment pairs.  
%Then, a graph-based kernel can be obtained with a two-phase approach: constructing distance matrix, and then kernelizing the distance matrix.  In the first phase, given a road network $G$, we define the adjacency matrix $\mathcal{A}$ as
%\begin{equation}
 % a_{ss'} =
  %\left\{\begin{array}{cl}
   % w(s,s') & \mbox{if}\ (s,s')\in E \vee (s',s)\in E, \\
%0 & \mbox{otherwise}.
%\end{array}\right.
%\label{eq:a}
%\end{equation}
%(Operator $[.]$ is used to manipulate a specific element in a matrix or vector.)
%The weight function models the dissimilarity of traffic phenomena between two adjacent segments. To quantify the dissimilarity, we can use
%We can view this dissimilarity as the `virtual' distance between adjacent segments.  
The next step is to compute the shortest path distance $d(s,s')$ between all pairs of road segments $s,s'\in V$ (i.e., using Floyd-Warshall or Johnson's algorithm) with respect to the topology of the weighted directed graph $G$.
% and store the results in a distance matrix denoted by $\mathcal{G}$.
Such a distance function is again a measure of dissimilarity, rather than one of similarity, as required by a kernel function. 
%To obtain the distance for any arbitrary segment pair, we can find all-pair shortest path based on $\mathcal{A}$. The weights of all-shortest paths are stored in a distance matrix denoted by $\mathcal{G}$.
% how to obtain kernel from distance?
%The second phase is to model $\Sigma(.,.;G)$ by kernelizing distance matrix $\mathcal{G}$.  
%As known, kernel function of GP models the correlation which is a similarity measure, while distance value in $\mathcal{G}$ is a measure of dissimilarity.  
Furthermore, a valid GP kernel needs to be positive semidefinite and symmetric \citep{scholkopf02}, which are clearly violated by $d$.  
%Clearly, $\mathcal{G}$ also does not satisfy this requirement.  

To construct a valid GP kernel from $d$, multi-dimensional scaling \citep{borg05} is
applied to embed the domain of road segments into the $p'$-dimensional Euclidean space $\mathbb{R}^{p'}$.
%$d$-dimensional Euclidean space $\mathcal{X}$.  
Specifically, a mapping $g: V\rightarrow\mathbb{R}^{p'}$ is determined by minimizing the squared loss
%
%\begin{equation*}
$  g^\ast=\mathop{\arg\min}_g \sum_{s,s'\in V}
  (d(s,s')-\|g(s)-g(s')\|)^2
$.
%  \label{mds}
%\end{equation*}
%
%Let $x\triangleq g^\ast (s)$ and $x'\triangleq g^\ast (s')$ for any $s,s'\in V$.
With a small squared loss, the Euclidean distance $\|g^{\ast}(s)-g^{\ast}(s')\|$ between $g^{\ast}(s)$ and $g^{\ast}(s')$ is expected to closely approximate the shortest path distance $d(s,s')$ between any pair of road segments $s$ and $s'$.
%, where $x_i,x_j\in \mathcal{X}$.  
After embedding into Euclidean space, a conventional kernel function such as the squared exponential one \citep{rasmussen06} can be used:\vspace{-2.5mm}
%
\begin{equation*}
  k(s,s')=\sigma_s^2\exp\left({-\frac{1}{2}\sum_{i=1}^{p'}\left(\frac{[g^{\ast}(s)]_i-[g^{\ast}(s')]_i}{\ell_i}\right)^2}\right)
\label{gkernel}\vspace{-2.5mm}
\end{equation*}
%
where $[g^{\ast}(s)]_i$ ($[g^{\ast}(s')]_i$) is the $i$-th component of the $p'$-dimensional vector $g^{\ast}(s)$ ($g^{\ast}(s')$), and the hyperparameters $\sigma_s$, $\ell_1, \ldots, \ell_{p'}$ are, respectively, signal variance and length-scales that can be learned using maximum likelihood estimation \citep{rasmussen06}. The resulting kernel function $k$\footnote{For spatiotemporal traffic modeling, the kernel function $k$ can be extended to account for the temporal dimension.} is guaranteed to be valid. \vspace{-2mm}
%$$\Sigma(\mathbf{s_i},\mathbf{s_j})$$

\subsection {Subset of Data Approximation}
\label{sec:sgp}\vspace{-3mm}
%Drawback of GP
Although the GP is an effective predictive model, it faces a practical limitation of cubic time complexity in the number $|D|$ of observations; this can be observed from computing the posterior distribution (i.e., (\ref{pmu}) and (\ref{pcov})), which requires inverting covariance matrix $\Sigma_{DD}$ and incurs $\set{O}(|D|^3)$ time.
%GP is an effective model for active data gathering with mobile sensors, however, to compute the predictive distribution, the kernel matrix  should be inversed in  given observed data $D$.
If $|D|$ is expected to be large, GP prediction cannot be performed in real time.  
%This is impractical for real world applications as $D$ could keep increasing to a large amount.  
For practical usage, we have to resort to computationally cheaper approximate GP prediction.
%Therefore, approximation of GP which is  is necessary for practical use.  

A simple method of approximation is to select only a subset $U$ of the entire set $D$ of observed road segments (i.e., $U\subset D$)
%, $|U|\ll |D|$) 
to compute the posterior distribution of the measurements at any set $Y\subseteq V\setminus D$ of unobserved road segments.
%, rather than the full history of data.
Such a \emph{subset of data} (SoD) approximation method produces the following predictive Gaussian distribution, which closely resembles that of the full GP model (i.e., by simply replacing $D$ in (\ref{pmu}) and (\ref{pcov}) with $U$):\vspace{-2mm}
% is called a sparse GP because the covariance matrix in a sparse GP is more sparse as compared to the original one.  
%Then, the predictive distribution of sparse GP using  on an unobserved set $Y$ is specified by
%
\begin{equation}
%  \mu_{Y}^{\mbox{\scriptsize SoD}}=\mu_Y+\Sigma_{YU}\Sigma_{UU}^{-1} (z_U - \mu_U) 
  \mu_{Y|U}=\mu_Y+\Sigma_{YU}\Sigma_{UU}^{-1} (z_U - \mu_U) 
\label{eq:fsod}\vspace{-2mm}
\end{equation}
\begin{equation}
%  \Sigma_{YY}^{\mbox{\scriptsize SoD}}=\Sigma_{YY}-\Sigma_{YU}\Sigma_{UU}^{-1}\Sigma_{UY}\ .
  \Sigma_{YY|U}=\Sigma_{YY}-\Sigma_{YU}\Sigma_{UU}^{-1}\Sigma_{UY}\ .
\label{pcovsod}\vspace{-1mm}
\end{equation}
%
Notice that the covariance matrix $\Sigma_{UU}$ to be inverted only incurs $\set{O}(|U|^3)$ time, which is independent of $|D|$.  
%The time complexity of the SoD approximation method reduces to $\set{O}(|U|^3)$, which is independent of $|D|$.  

The predictive performance of SoD approximation is sensitive to the selection of subset $U$.
%depends on which data points are selected into $U$.  Intuitively, 
In practice, random subset selection often yields poor performance.  
This issue can be resolved by actively selecting an informative subset $U$ in an iterative greedy manner:
%actively selecting the data points iteratively.  A greedy approach is briefly described as follows:  
Firstly, $U$ is initialized to be an empty set. 
Then, all road segments in $D\setminus U$ are scored based on a criterion that can be chosen from, for example, the works of \citep{Guestrin08,lawrence03,seeger03}.
% a specific criterion to be discussed later. 
The highest-scored segment is selected for inclusion in $U$ and removed from $D$.
%the selection procedure starts by computing the scores of data points in . According to the score, the best data point is selected into $U$ and removed from $D$.  
This greedy selection procedure is iterated until $U$ reaches a pre-defined size.  
%Note that the criterion can be any one introduced in the works.  
Among the various criteria introduced earlier, the differential entropy score \citep{lawrence03} is reported to perform well \citep{oh10}; it is a monotonic function of the posterior variance $\Sigma_{ss|U}$ (\ref{pcovsod}), thus resulting in the greedy selection of a segment $s\in D\setminus U$ with the largest variance in each iteration. \vspace{-3mm}
%$\mathbb{H}[Z_s|Z_{D\setminus U}]$
%Due to the greedy fashion, the entropy of a Gaussian random variable  with variance $\sigma^2$ is $\log(2\pi e\sigma^2)/2$ (\ref{entropy}), thus the IVM simply reduces to selecting the segment with the largest predictive variance in each iteration.
%
\section{Decentralized Data Fusion}
\label{sec:ddf}\vspace{-3mm}
%
In the previous section, two centralized data fusion approaches to exact (i.e., (\ref{pmu}) and (\ref{pcov})) and approximate (i.e., (\ref{eq:fsod}) and (\ref{pcovsod})) GP prediction are introduced.
% are introduced, namely, the full GP  in Section~\ref{sec:gpgk}) and the sparse SoD approximation methods  in Section~\ref{sec:sgp}).
%As we known, the FGP and SoD can be used to solve the ASMS problem, however, both methods have to work in centralized fashion. 
%To compute the joint entropy in (\ref{maxent}), both methods are exponential in the number of mobile sensors and length of planning horizon, which is undesirable for online prediction. 
In this section, we will discuss the decentralized data fusion component of our D$^2$FAS algorithm, which distributes the computational load among the mobile sensors to achieve efficient and scalable approximate GP prediction.
%Our approach relies on a different technique which still exploits the full data $D$ but uses a sparse summarized structure to preserve the efficiency.  

The intuition of our decentralized data fusion algorithm is as follows: Each of the $K$ mobile sensors constructs a local summary of the observations taken along its own path in the road network and communicates its local summary to every other sensor. Then, it assimilates the local summaries received from the other sensors into a globally consistent summary, which is exploited for predicting the traffic phenomenon as well as active sensing.
This intuition will be formally realized and described in the paragraphs below.

While exploring the road network, each mobile sensor summarizes its local observations taken along its path based on a common support set $U\subset V$ known to all the other sensors. Its local summary is defined as follows:\vspace{-1.5mm}
% we first define the local summary for a mobile sensor to summarize its local data given a support set $U$ (the same notation as the subset in SoD is used to facilitate comparison, however they are conceptually different). 
%
\begin{defi}[Local Summary]
Given a common support set $U\subset V$ known to all $K$ mobile sensors, a set $D_k \subset V$ of observed road segments and a column vector $z_{D_k}$  of corresponding measurements local to mobile sensor $k$, its local summary is defined as a tuple $(\dot{z}_U^k, \dot{\Sigma}_{UU}^k)$ where\vspace{-0.5mm}
%
\begin{equation}
\dot{z}_{U}^k \triangleq\Sigma_{UD_k} \Sigma^{-1}_{D_k D_k|U} (z_{D_k} - \mu_{D_k})
\label{lsf}\vspace{-0.5mm}
\end{equation}
\begin{equation}
\dot{\Sigma}_{UU}^k \triangleq\Sigma_{UD_k}\Sigma_{D_k D_k|U}^{-1}\Sigma_{D_k U}
\label{lsc}\vspace{-0.5mm}
\end{equation}
such that $\Sigma_{D_k D_k|U}$ is defined in a similar manner to $($\ref{pcovsod}$)$.
\label{def1}\vspace{-1.5mm}
\end{defi}
%
\emph{Remark}. Unlike SoD (Section~\ref{sec:sgp}), the support set $U$ of road segments does not have to be observed, since the local summary (i.e., (\ref{lsf}) and (\ref{lsc})) is independent of the corresponding measurements $z_U$. 
%are not required in computing the local summary (i.e., (\ref{lsf}) and (\ref{lsc})). 
So, $U$ does not need to be a subset of $D=\bigcup^K_{k=1}D_k$. To select an informative support set $U$ from the set $V$ of all possible segments in the road network, an offline active selection procedure similar to that in the last paragraph of Section~\ref{sec:sgp} can be performed just once prior to observing data to determine $U$. In contrast, SoD has to perform online active selection every time new road segments are being observed.
%Notice that, unlike SoD, the measurements of $U$ in DADGF are not required to compute local summary, thus $U$ is not necessarily a subset of $D$.  This property enables the support set to be actively selected in an offline fashion.  To lower the performance reduction due to dropped data in $D/U$, SoD employs active selection to find a better representative subset. However, the active selection is computationally intensive and has to be performed online. 
%Since there are no measurements corresponding with $U$, the matrix $\Sigma_{D_i,D_i}|U$ is computed with zero noise by (\ref{pcov}).
%

By communicating its local summary to every other sensor, each mobile sensor can then construct a globally consistent summary from the received local summaries:\vspace{-1.5mm}
%Like the full data $D$ in FGP or SoD, the DADGF also needs to maintain a global view among mobile sensors to make consistent prediction. We call this view as global summary which is constructed by assimilating the local summaries distributed among $K$ mobile sensors. 
\begin{defi}[Global Summary] Given a common support set $U\subset V$ known to all $K$ mobile sensors and the local summary $(\dot{z}_U^k, \dot{\Sigma}_{UU}^k)$ of every mobile sensor $k = 1,\ldots,K$,
the global summary 
  %denoted as $\bar{S}_U$ 
is defined as a tuple $(\ddot{z}_U, \ddot{\Sigma}_{UU})$ where\vspace{-2mm}
%
\begin{equation}
\ddot{z}_U \triangleq\sum_{k=1}^{K}\dot{z}_{U}^k
\label{gsf}\vspace{-2mm}
\end{equation}
\begin{equation}
\ddot{\Sigma}_{UU} \triangleq\Sigma_{UU}+\sum_{k=1}^{K} \dot{\Sigma}_{UU}^k \ .
\label{gsc}\vspace{-6mm}
\end{equation}
\label{def2}
\end{defi}
\emph{Remark}. In this paper, we assume all-to-all communication between the $K$ mobile sensors.
% with $O(K^2)$ communication complexity which is the same for all three approaches.  
Supposing this is not possible and each sensor can only communicate locally with its neighbors, the summation structure of the global summary (specifically, (\ref{gsf}) and (\ref{gsc})) makes it amenable to be constructed using distributed consensus filters \citep{saber05a}. We omit these details since they are beyond the scope of this paper.

Finally, the global summary is exploited by each mobile sensor to compute a globally consistent predictive Gaussian distribution, as detailed in Theorem~\ref{thmdad}A below, as well as to perform decentralized active sensing (Section~\ref{sec:dadgf}):\vspace{-2mm}
%
\begin{thm}
Let a common support set $U\subset V$ be known to all $K$ mobile sensors.\vspace{-4mm}
%  With $K$ mobile sensors, if the global summary is obtained by (\ref{gsf}) and (\ref{gsc}) then for an arbitrary mobile sensor
\begin{description}
\item[A.\ \ ] Given the global summary $(\ddot{z}_U, \ddot{\Sigma}_{UU})$, each mobile sensor computes a globally consistent predictive Gaussian distribution
  $\mathcal{N}(\overline{\mu}_{Y},\overline{\Sigma}_{YY})$ of the measurements at any set $Y$ of unobserved road segments where\vspace{-1mm}
% can be found locally such that
%
\begin{equation}
    \overline{\mu}_{Y} \triangleq\mu_Y + \Sigma_{YU}\ddot{\Sigma}_{UU}^{-1}\ddot{z}_U
\label{mupitc}\vspace{-1mm}
\end{equation}
\begin{equation}
  \overline{\Sigma}_{YY} \triangleq\Sigma_{YY}-\Sigma_{YU}(\Sigma_{UU}^{-1}-\ddot{\Sigma}
_{UU}^{-1})\Sigma_{UY} \ .
\label{varpitc}\vspace{-0mm}
\end{equation}
% 
\item[B.\ \ ] Let $\mathcal{N}(\mu_{Y|D}^{\mbox{\tiny\emph{PITC}}},\Sigma_{YY|D}^{\mbox{\tiny\emph{PITC}}})$ be  the predictive Gaussian distribution computed by the centralized sparse partially independent training conditional (PITC) approximation of GP model \citep{candela05} where\vspace{-2mm}
%with support set $U$ and data $D$ partitioned into $K$ blocks.
\begin{equation} 
  \mu_{Y|D}^{\mbox{\tiny\emph{PITC}}} \triangleq\mu_Y+\Gamma_{YD}\left(\Gamma_{DD}+\Lambda\right)^{-1}\hspace{-1mm}\left(z_D-\mu_D \right)
  \label{eq:pitc.mu}\vspace{-2mm}
\end{equation}   
\begin{equation}   
\Sigma_{YY|D}^{\mbox{\tiny\emph{PITC}}} \triangleq\Sigma_{YY}-\Gamma_{YD}\left(\Gamma_{DD}+\Lambda\right)^{-1}\Gamma_{DY}
  \label{eq:pitc.var} \vspace{-0mm}
\end{equation} 
such that \vspace{-0mm}
\begin{equation} 
  \Gamma_{BB'} \triangleq\Sigma_{BU}\Sigma_{UU}^{-1}\Sigma_{UB'} 
  \label{eq:sparseq}
\end{equation} 
and 
%$U$ is the support set for sparse GP models.
$\Lambda$
%\triangleq
%\mbox{blockdiag}[\Sigma_{DD|U}]=
%\mbox{blockdiag}[\Sigma_{DD}-\Gamma_{DD}]$ 
is a 
%$|D|\times |D|$ 
block-diagonal matrix constructed from the $K$ diagonal blocks of $\Sigma_{DD|U}$,
%=\Sigma_{DD}-\Gamma_{DD}$, 
each of which is a 
%$|D_k|\times |D_k|$ 
matrix $\Sigma_{D_k D_k|U}$
%=\Sigma_{D_k D_k}-Q_{D_k D_k}$ 
for $k=1,\ldots,K$ where $D=\bigcup_{k=1}^K D_k$.
%\Lambda_k \triangleq
Then, $\overline{\mu}_{Y}\hspace{-0.5mm}=\hspace{-0.5mm}\mu_{Y|D}^{\mbox{\tiny\emph{PITC}}}$ and $\overline{\Sigma}_{YY}\hspace{-0.5mm}=\hspace{-0.5mm}\Sigma_{YY|D}^{\mbox{\tiny\emph{PITC}}}$.
\end{description} 
%
\label{thmdad}\vspace{-4mm}
\end{thm}
%
% (i) is valid immediately, because the segments in $U$ and $Y$ are given locally. 
\ifthenelse{\value{sol}=1}{The proof of Theorem~\ref{thmdad}B is given in Appendix~\ref{equ}.} 
{Its proof is given in \citep{LowArxiv12}.} 
The equivalence result of Theorem~\ref{thmdad}B bears two implications:

\emph{Remark} $1$.
The computation of PITC can be parallelized and distributed among the mobile sensors in a Google-like MapReduce paradigm \citep{Chu07b}, thereby improving the time efficiency of prediction: Each of the $K$ mappers (sensors) is tasked to compute its local summary while the reducer (any sensor) sums these local summaries into a global summary, which is then used to compute the predictive Gaussian distribution.
%The computational load of the centralized PITC approximation of GP model can be distributed among $K$ mobile sensors, thereby improving the time efficiency of prediction. 
Supposing $|Y|\leq |U|$ for simplicity, the $\set{O}\hspace{-1mm}\left(|D|((|D|/K)^2+|U|^2)\right)$ time incurred by PITC can be reduced to $\set{O}\hspace{-1mm}\left((|D|/K)^3+|U|^3+|U|^2 K\right)$ time of running our decentralized algorithm on each of the $K$ sensors, the latter of which scales better with increasing number $|D|$ of observations.

\emph{Remark} $2$.
%The implication of Theorem~\ref{thmdad}B is that since the predictive Gaussian distributions of D$^2$FAS and PITC are equivalent, 
We can draw insights from PITC to elucidate an underlying property of our decentralized algorithm: 
It is assumed that $Z_{D_1},\ldots,Z_{D_K}, Z_Y$ are conditionally independent given the measurements at the support set $U$ of road segments.
To potentially reduce the degree of violation of this assumption, an informative support set $U$ is actively selected, as described earlier in this section.
Furthermore, the experimental results on real-world urban road network data\footnote{\cite{candela05} only illustrated the predictive performance of PITC on a simulated toy example.} (Section~\ref{sec:exp})
show that D$^2$FAS can achieve predictive performance comparable to that of the full GP model while enjoying significantly lower computational cost, thus demonstrating the practicality of such an assumption for predicting traffic phenomena.
The predictive performance of D$^2$FAS can be improved by increasing the size of $U$ at the expense of greater time and communication overhead. \vspace{-2mm}
%trade off time and communication for better predictive power with larger $U$
%$P(Z_D,Z_Y|z_u) = P(Z_D|z_u) P(Z_Y|z_u)$
%$P(Z_D|z_u) = P(Z_{D_1},\ldots,Z_{D_K}|z_u) = \prod^K_{k=1} P(Z_{D_k}|z_u)$
% Sparse posterior covariance matrix given the support set?
%
\section{Decentralized Active Sensing}
\label{sec:dadgf}\vspace{-3mm}
%overview
%In this section, we first formulate the problem of active sensing with mobile sensors and show that FGP and SoD can solve the problem in centralized fashion. Then, we present our novel decentralized active data fusion and gathering (DADGF) algorithm.  
%\subsection{Active Sensing with Mobile Sensors}
The problem of active sensing with $K$ mobile sensors is formulated as follows: 
Given the set $D_k\subset V$ of observed road segments and the currently traversed road segment $s_k\in V$ of every mobile sensor $k =1,\ldots,K$, the mobile sensors have to coordinate to select the most informative walks $w^{\ast}_1,\ldots, w^{\ast}_K$ of length (i.e., number of road segments) $L$ each and with respective origins $s_1,\ldots,s_K$ in the road network $G$: \vspace{-1mm}
% Assuming we have $K$ mobile sensors, any arbitrary sensor $k$ maintains knowledge such as observed segments $D_k\subset V$ and current position $p_k\in V$.  At certain time point, each mobile sensor is planned to explore a trajectory of length $L$. The set of candidate trajectories is denoted by $C_k$, where $\mathbf{c_k}\in C_k$ is an $L$-dimensional vector.  Clearly, $C_k$ is constrained by $p_k$ and topology of $G$. The mobile sensors evaluate the joint conditional entropy to determine the most informative paths $(\mathbf{c_1^*},\dots,\mathbf{c_K^*})$ by 
%
\begin{equation}
\hspace{-0.3mm}
\displaystyle(w^{\ast}_1,\ldots, w^{\ast}_K)\hspace{-0.8mm} = \hspace{-0.8mm}\mathop{\arg\max}_{(w_1,\ldots, w_K)}\mathbb{H}\hspace{-1mm}\left[Z_{\bigcup^K_{k=1} Y_{w_k}} \Big{|} Z_{\bigcup^K_{k=1} D_k}\right]\vspace{-1mm}
% =   \displaystyle\mathop{\arg\max}_{(w_1,\ldots, w_K)} \left|\Sigma_{\bigcup^K_{k=1} Y_{w_k} \bigcup^K_{k=1} Y_{w_k}|\bigcup^K_{k=1} D_k}\right|
\label{maxent}
\end{equation}
%
where $Y_{w_k}$ denotes the set of unobserved road segments induced by the walk $w_k$.
To simplify notation, let a joint walk be denoted by $w\triangleq (w_1,\ldots,w_K)$ (similarly, for $w^\ast$) and its induced set of unobserved road segments be $Y_{w}\triangleq \bigcup^K_{k=1} Y_{w_k}$ from now on.
Interestingly, it can be shown using the chain rule for entropy that these maximum-entropy walks $w^{\ast}$ minimize the posterior joint entropy (i.e., $\mathbb{H}[Z_{V\setminus(D\bigcup Y_{w^{\ast}})} | Z_{D\bigcup Y_{w^{\ast}}}]$) of the measurements at the remaining unobserved segments (i.e., $V\setminus(D\bigcup Y_{w^{\ast}})$) in the road network.
% where $D=\bigcup_{k=1}^K D_k$.
After executing the walk $w^{\ast}_k$, each mobile sensor $k$ observes the set $Y_{w^{\ast}_k}$ of road segments and updates its local information:\vspace{-0.8mm}
%Then mobile sensor $k$ gathers traffic data on $Y_k$ which contains the segments along $c_k^*$, and updates the local knowledge:
%
\begin{equation}
  D_k\leftarrow D_k\bigcup Y_{w^{\ast}_k} \ , z_{D_k}\leftarrow z_{D_k\bigcup Y_{w^{\ast}_k}}, s_k\leftarrow \mbox{terminus of}\ w^{\ast}_k\ .\vspace{-0mm}
\label{upDk}
\end{equation}
%\begin{equation}
%
%\label{upPk}
%\end{equation}
%
%If necessary, the active sensing problem (\ref{maxent}) can be solved repeatedly to produce subsequent maximum-entropy walks.
Evaluating the Gaussian posterior entropy term in (\ref{maxent}) involves computing a posterior covariance matrix (\ref{entropy}) using one of the data fusion methods described earlier:
If (\ref{pcov}) of full GP model (Section~\ref{sec:gpgk}) or (\ref{pcovsod}) of SoD (Section~\ref{sec:sgp}) is to be used, then the observations that are gathered distributedly by the sensors have to be fully communicated to a central data fusion center.
In contrast, our decentralized data fusion algorithm (Section~\ref{sec:ddf}) only requires communicating local summaries (Definition~\ref{def1}) to compute (\ref{varpitc}) for solving the active sensing problem (\ref{maxent}):\vspace{-0.5mm}
\begin{equation}
  \displaystyle w^{\ast}=\displaystyle\mathop{\arg\max}_{w}\ \overline{\mathbb{H}}\hspace{-0.5mm}\left[Z_{Y_w}\right]\ ,\label{maxent2}\vspace{-1.5mm}
\end{equation}
\begin{equation}
\overline{\mathbb{H}}\hspace{-0.5mm}\left[Z_{Y_w}\right]\triangleq\frac{1}{2}\log\hspace{-0.5mm}\left(2\pi e\right)^{\left|Y_w\right|}\left|\overline{\Sigma}_{Y_{w} Y_{w}}\right|\ .\label{maxent21}\vspace{-1mm}
\end{equation}
%If full GP or SoD is to be performed separately on each of the $K$ mobile sensors rather than centrally, then the observations that are gathered distributedly by the sensors have to be fully communicated to every sensor.
%In contrast, D$^2$FAS only requires exchanging local summaries (Definition~\ref{def1}) between sensors.
%consolidated centrally
%, after which the posterior entropy term in (\ref{maxent}) can be computed using (\ref{pcov}) \& (\ref{entropy}) and  (\ref{pcovsod}) \& (\ref{entropy}), respectively.
%From (\ref{dmaxent}), we can see that the exploring strategy can be determined in decentralized fashion. 
Without imposing any structural assumption, solving the active sensing problem (\ref{maxent2}) will be  prohibitively expensive due to the space of possible joint walks $w$ that grows exponentially in the number $K$ of mobile sensors.
%, thus severely compromising the scalability. 
%As a result, solving the active sensing problem (\ref{maxent}) with FGP or SoD incurs exponential time in the number $K$ of mobile sensors, thus severely compromising their scalability. 
%Their time and communication complexities are shown in Table~\ref{comp}.
%The full GP (FGP) model (Section \ref{sec:gpgk}) can be used with mobile sensors to solve ASMS problem.  As FGP depends on the full data which are distributed among $K$ mobile sensors, a global consistent set $D$ is first constructed from local data by $D=\bigcup_{k=1}^K D_k$.
%Similarly, SoD model (Section~\ref{sec:sgp}) is also applicable with the differences that a subset $U$ is selected actively from $D$ first, and then (\ref{maxent}) is evaluated by (\ref{pcovsod}) and (\ref{entropy}) instead. 
%
%\subsection{Decentralized Active Data Gathering and Fusion}
%\label{sec:dadgf}
%As we known, the FGP and SoD can be used to solve the ASMS problem, however, both methods have to work in centralized fashion. To compute the joint entropy in (\ref{maxent}), both methods are exponential in the number of mobile sensors and length of planning horizon, which is undesirable for online prediction. In this section, we present our novel decentralized approach called Decentralized Active Data Gathering and Fusion (DADGF).   Finally, we notice that, to find the maximum entropy strategy by (\ref{maxent}), the number of joint entropy evaluation grows exponentially in the number of mobile sensors $K$. 
%To overcome the same scalability issue with solving the active sensing problem (\ref{maxent}) using D$^2$FAS, we extend an earlier assumption in Section~\ref{sec:ddf} that $Z_{D_1},\ldots,Z_{D_K}, Z_{Y_{w_1}},\ldots, Z_{Y_{w_K}}$ are conditionally independent given the measurements at the support set $U$ of road segments.
%Therefore, we make assumption that $P(\tilde{f}_Y|\overline{S}_U)=\prod_{i=1}^K P(\tilde{f}_{Y_i}|\overline{S}_U)$, where $Y_i$ consists of the segments in $c_i \in C_i$ and $Y=\cup_{i=1}^K Y_i$.  Then
To overcome this scalability issue for D$^2$FAS, our key idea is to construct a block-diagonal matrix whose log-determinant closely approximates that of $\overline{\Sigma}_{Y_{w} Y_{w}}$ (\ref{varpitc}) and exploit the property that the log-determinant of such a block-diagonal matrix can be decomposed into a sum of log-determinants of its diagonal blocks, each of which depends only on the walks of a disjoint subset of the $K$ mobile sensors.
Consequently, the active sensing problem can be partially decentralized, leading to a reduced space of possible joint walks to be searched, as detailed in the rest of this section. 

Firstly, we extend an earlier assumption in Section~\ref{sec:ddf}:
$Z_{D_1},\ldots,Z_{D_K}, Z_{Y_{w_1}},\ldots,Z_{Y_{w_K}}$ are conditionally independent given the measurements at the support set $U$ of road segments. 
Then, it can be shown via the equivalence to PITC (Theorem~\ref{thmdad}B) 
that $\overline{\Sigma}_{Y_{w} Y_{w}}$ (\ref{varpitc}) comprises diagonal blocks of the form $\overline{\Sigma}_{Y_{w_k} Y_{w_k}}$ for $k=1,\ldots,K$ and off-diagonal blocks of the form $\Sigma_{Y_{w_k}U}\ddot{\Sigma}
_{UU}^{-1}\Sigma_{UY_{w_{k'}}}$ for $k,k'=1,\ldots,K$ and $k\neq k'$.
In particular, each off-diagonal block of $\overline{\Sigma}_{Y_{w} Y_{w}}$ represents the correlation of measurements between the unobserved road segments $Y_{w_k}$ and $Y_{w_{k'}}$ along the respective walks $w_k$ of sensor $k$ and $w_{k'}$ of sensor $k'$. If the correlation between some pair of their possible walks is high enough, then
their walks have to be coordinated.
%A high enough correlation should  inspire the coordination of the walks of sensors $k$ and $k'$.
%Using the components in the off-diagonal blocks of $\overline{\Sigma}_{Y_{w} Y_{w}}$, 
This is formally realized by the following coordination graph over the $K$ sensors:\vspace{-1.5mm}
\begin{defi}[Coordination Graph]\label{coord}
Define the coordination graph to be an undirected graph $\mathcal{G}\triangleq (\mathcal{V}, \mathcal{E})$ that comprises\vspace{-2mm}
\squishlisttwo
\item a set $\mathcal{V}$ of vertices denoting the $K$ mobile sensors, and 
\item a set $\mathcal{E}$ of edges denoting coordination dependencies between sensors such that there exists an edge $\{k,k'\}$ incident with sensors $k\in \mathcal{V}$ and $k'\in \mathcal{V}\setminus\{k\}$ iff\vspace{-1mm}
\begin{equation}
\max_{s\in Y_{W_k}, s'\in Y_{W_{k'}}} \left|\Sigma_{sU}\ddot{\Sigma}_{UU}^{-1}\Sigma_{Us'}\right| > \varepsilon
\label{rel}\vspace{-1mm}
\end{equation}
for a predefined constant $\varepsilon >0$ where $W_k$ denotes the set of possible walks of length $L$ of mobile sensor $k$ from origin $s_k$ in the road network $G$ and $Y_{W_k} \triangleq \bigcup_{w_k\in W_k}Y_{w_k}$.\vspace{-1mm}
\squishend
\end{defi}
%Intuitively, each edge indicates coordination between sensors $k$ and $k'$ when some pair of their walks contain highly correlated road segments.
\emph{Remark}. The construction of $\mathcal{G}$ can be decentralized as follows:
Since $\ddot{\Sigma}_{UU}$ is symmetric and positive definite,
it can be decomposed by Cholesky factorization into $\ddot{\Sigma}_{UU} = \Psi \Psi^{\top}$ where $\Psi$ is a lower triangular matrix and $\Psi^{\top}$ is the transpose of $\Psi$.
Then,
%each component in any off-diagonal block of $\overline{\Sigma}_{Y_{w} Y_{w}}$ can be factorized into 
$\Sigma_{sU}\ddot{\Sigma}_{UU}^{-1}\Sigma_{Us'} =  (\Psi\backslash\Sigma_{Us})^{\top}\Psi\backslash\Sigma_{Us'}$ where $\Psi\backslash B$ denotes the column vector $\phi$ solving $\Psi\phi = B$.
That is, $\Sigma_{sU}\ddot{\Sigma}_{UU}^{-1}\Sigma_{Us'}$ (\ref{rel}) can be expressed as a dot product of two vectors $\Psi\backslash\Sigma_{Us}$ and $\Psi\backslash\Sigma_{Us'}$; this property is exploited to determine adjacency between sensors in a decentralized manner:\vspace{-2mm} 
\begin{defi}[Adjacency]\label{adj}
Let\vspace{-1.5mm}
\begin{equation}
\Phi_k \triangleq\{\Psi\backslash\Sigma_{Us}\}_{s\in Y_{W_k}} 
\label{phik}\vspace{-1.5mm}
\end{equation}
for $k=1,\ldots,K$.
A sensor $k\in \mathcal{V}$ is adjacent to sensor $k'\in \mathcal{V}\setminus\{k\}$ in coordination graph $\mathcal{G}$ iff\vspace{-1.5mm}
\begin{equation}
\max_{\phi\in \Phi_k, \phi'\in \Phi_{k'}} \left|\phi^{\top} \phi'\right| > \varepsilon \ .
\label{rel2}\vspace{-3.0mm}
\end{equation}
\end{defi}
It follows from the above definition that if each sensor $k$ constructs $\Phi_k$ and exchanges it with every other sensor, then it can determine its adjacency to all the other sensors and store this information in a column vector $a_k$ of length $K$ with its $k'$-th component being defined as follows:\vspace{-2mm}
\begin{equation}
\hspace{-1mm}
\left[a_k\right]_{k'}\hspace{-0.1mm} =\hspace{-0.1mm}\Bigg{\{}
\begin{array}{rl}
1 & \mbox{if sensor $k$ is adjacent to sensor $k'$,}\vspace{1mm}\\
0 & \mbox{otherwise.}
\end{array}
\label{adjv}\vspace{-2mm}
\end{equation}
%its $k'$-th component is of value $1$ if sensor $k$ is adjacent to sensor $k'$, and $0$ otherwise.
By exchanging its adjacency vector $a_k$ with every other sensor, each sensor can construct a globally consistent adjacency matrix $A_{\mathcal{G}} \triangleq (a_1\ldots a_K)$ to represent coordination graph $\mathcal{G}$.
%After receiving $\Phi_{k'}$ for $k'\in \mathcal{V}\setminus \{k\}$, we know using (\ref{rel}) that sensor $k$ is adjacent to sensor $k'$ iff
%Each sensor can compute its own adjacency information and exchange it with every other sensor to construct a globally consistent adjacency matrix $A_{\mathcal{G}}$.
%Partition $\mathcal{V}$ into a collection of $\mathcal{K}$ disjoint subsets $\{ \mathcal{V}_i \}_{i=1,\ldots,\mathcal{K}}$.

Next, by computing the connected components (say, $\mathcal{K}$ of them) of coordination graph $\mathcal{G}$, their resulting vertex sets partition the set $\mathcal{V}$ of $K$ sensors into $\mathcal{K}$ disjoint subsets $\mathcal{V}_1,\ldots,\mathcal{V}_{\mathcal{K}}$ such that the sensors within each subset have to coordinate their walks.
Each sensor can determine its residing connected component in a decentralized way by performing a depth-first search in $\mathcal{G}$ starting from it as root.
%each sensor $k$ can run depth-first search with root $k$ to determine the connected component it resides in.

Finally, construct a block-diagonal matrix $\widehat{\Sigma}_{Y_{w} Y_{w}}$ to comprise diagonal blocks of the form $\overline{\Sigma}_{Y_{w_{\mathcal{V}_n}} Y_{w_{\mathcal{V}_n}}}$ for $n = 1,\ldots,\mathcal{K}$
where $w_{\mathcal{V}_n}\triangleq (w_k)_{k\in \mathcal{V}_n}$ and $Y_{w_{\mathcal{V}_n}}\triangleq\bigcup_{k\in \mathcal{V}_n} Y_{w_k}$.
%Such an assumption is not uncommon: it is often made in order to calculate the widely-used sum of posterior variances (i.e., mean-squared error) criterion (Section~\ref{sec:gpgk}).
%In practice, this assumption usually becomes less restrictive when the number $|\set{D}|$ of observed road segments increases to potentially reduce the degree of violation of conditional independence, the correlation of measurements between road segments decreases, and/or the mobile sensors are sufficiently far apart. 
%Using the chain rule for entropy and subsequently the conditional independence assumption, the active sensing problem (\ref{maxent}) reduces to 
The active sensing problem (\ref{maxent2}) is then approximated by\vspace{-2mm}
\begin{equation}
\hspace{-0.5mm}
\begin{array}{l}
\displaystyle\max_{w}\ \frac{1}{2}\log\hspace{-0.5mm}\left(2\pi e\right)^{\left|Y_w\right|}\left|\widehat{\Sigma}_{Y_{w} Y_{w}}\right|\\ 
\equiv\displaystyle\max_{(w_{\mathcal{V}_1},\ldots, w_{\mathcal{V}_{\mathcal{K}}})}\sum_{n=1}^{\mathcal{K}}\log\hspace{-0.5mm}\left(2\pi e\right)^{\left|Y_{w_{\mathcal{V}_n}}\right|}\left|\overline{\Sigma}_{Y_{w_{\mathcal{V}_n}} Y_{w_{\mathcal{V}_n}}}\right|\\
=\displaystyle\sum_{n=1}^{\mathcal{K}}\max_{w_{\mathcal{V}_n}}\ \log\hspace{-0.5mm}\left(2\pi e\right)^{\left|Y_{w_{\mathcal{V}_n}}\right|}\left|\overline{\Sigma}_{Y_{w_{\mathcal{V}_n}} Y_{w_{\mathcal{V}_n}}}\right|,
\end{array}
\label{bdiaga}\vspace{-2mm}
\end{equation}
which can be solved in a partially decentralized manner by each disjoint subset $\mathcal{V}_n$ of mobile sensors:\vspace{-2mm}
%a decentralized active sensing problem results: 
%
\begin{equation}
%\begin{array}{c}
%\hspace{-3mm}
\displaystyle \widehat{w}_{\mathcal{V}_n} = 
\mathop{\arg\max}_{w_{\mathcal{V}_n}}\ \log\hspace{-0.5mm}\left(2\pi e\right)^{\left|Y_{w_{\mathcal{V}_n}}\right|}\left|\overline{\Sigma}_{Y_{w_{\mathcal{V}_n}} Y_{w_{\mathcal{V}_n}}}\right|.
%\\
 % \displaystyle\widehat{\mathbb{H}}[Z_{Y_{w_k}}]\triangleq \frac{1}{2}\log(2\pi e)^{\left|Y_{w_k}\right|}\left|\widehat{\Sigma}_{Y_{w_k}Y_{w_k}}\right|
 % \end{array}
  \label{dmaxent}\vspace{-2mm}
\end{equation}
Our active sensing algorithm becomes fully decentralized if $\varepsilon$ is set to be sufficiently large: more sensors become isolated in $\mathcal{G}$, consequently decreasing the size $\displaystyle\kappa\triangleq\max_{n}\left|\mathcal{V}_n\right|$ of its largest connected component to $1$.
As shown in Section~\ref{sec:time}, 
decreasing $\kappa$ improves its time efficiency.
On the other hand, it tends to a centralized behavior (\ref{maxent2}) by setting $\varepsilon\rightarrow 0^{+}$: $\mathcal{G}$ becomes near-complete, thus resulting in $\kappa\rightarrow K$. \vspace{-2mm}
% (Definition~\ref{coord})
% and $\mathcal{K}\rightarrow 1$.
%In practice, the degree of decentralization can be controlled by varying $\varepsilon$:
%as $\varepsilon\rightarrow 0^{+}$,  Consequently, our active sensing algorithm tends to a .
%On the other hand, as $\varepsilon\rightarrow\infty$, $\kappa\rightarrow 1$ and $\mathcal{K}\rightarrow K$. As a result, our algorithm tends to a fully decentralized behavior.
%It turns out that these same conditions control the time efficiency.
% In practice, (a) increasing K or L tends to increase kappa with fixed epsilon; (b) decreasing epsilon tends to increase kappa to K and decrease \mathcal{K} to 1, centralized; (c) increasing epsilon tends to decrease kappa to 1 and increase \mathcal{K} to K, fully decentralized
%
%To ease notations, let $\widehat{w}$.
%such that the second equality follows from (\ref{entropy}) and the posterior covariance matrix $\Sigma_{Y_{w_k}Y_{w_k}|D}$ can be obtained using one of the data fusion methods described earlier, specifically, using (\ref{pcov}) of full GP model (Section~\ref{sec:gpgk}), (\ref{pcovsod}) of SoD (Section~\ref{sec:sgp}), or (\ref{varpitc}) of D$^2$FAS (Section~\ref{sec:ddf}).
%
%where $\mathbb{H}[Z_{Y_{w_k}} | Z_{D}]$ is defined in a similar way to (\ref{entropy}).
%it involves evaluating Gaussian posterior entropies, which in turn require computing predictive Gaussian distributions (in particular, covariance matrices) using one of the data fusion methods in Sections~\ref{sec:gpgk} and~\ref{sec:ddf}. 

Let
\vspace{-3.5mm} 
\begin{equation}
\displaystyle\xi\triangleq\max_{n,w_{\mathcal{V}_n},i,i'}\left|\left[\left(\overline{\Sigma}_{Y_{w_{\mathcal{V}_n}} Y_{w_{\mathcal{V}_n}}}\right)^{\hspace{-1mm}-1} \right]_{ii'}\right|
\label{max1}\vspace{-1mm}
\end{equation}
and $
\epsilon\triangleq\displaystyle 0.5\log 1\hspace{-0.5mm}\Big{/}\hspace{-1.5mm}\left(1\hspace{-0.5mm}-\hspace{-0.5mm}\left(K^{1.5}L^{2.5}\kappa\xi\varepsilon\right)^2\right)
$.
In the result below, we prove that the joint walk $\widehat{w}\triangleq(\widehat{w}_{\mathcal{V}_1},\ldots, \widehat{w}_{\mathcal{V}_{\mathcal{K}}})$ %induced by solving the partially decentralized active sensing problem (\ref{dmaxent}) 
is guaranteed to achieve an entropy $\overline{\mathbb{H}}\hspace{-0.5mm}\left[Z_{Y_{\widehat{w}}}\right]$ (i.e., by plugging $\widehat{w}$ into (\ref{maxent21})) that is not more than $\epsilon$ 
from the maximum entropy $\overline{\mathbb{H}}\hspace{-0.5mm}\left[Z_{Y_{w^{\ast}}}\right]$ achieved by joint walk $w^{\ast}$ (\ref{maxent2}):\vspace{-1mm}
\begin{thm}[Performance Guarantee]
If $K^{1.5}L^{2.5}\kappa\xi\varepsilon\hspace{-0.06mm}<1$, then
$\overline{\mathbb{H}}\hspace{-0.5mm}\left[Z_{Y_{w^{\ast}}}\right] - \overline{\mathbb{H}}\hspace{-0.5mm}\left[Z_{Y_{\widehat{w}}}\right]\leq\epsilon$.
%for $n=1,\ldots,\mathcal{K}$.
\label{thm2}\vspace{-1mm}
\end{thm}
\ifthenelse{\value{sol}=1}{The proof of Theorem~\ref{thm2} is given in Appendix~\ref{equ2}.} 
{Its proof is given in \citep{LowArxiv12}.} 
The implication of Theorem~\ref{thm2} is that our partially decentralized active sensing algorithm can perform comparatively well (i.e., small $\epsilon$) under the following favorable environmental conditions: (a) the network of $K$ sensors is not large, (b) length $L$ of each sensor's walk to be optimized is not long, (c) the largest subset of $\kappa$ sensors being formed to coordinate their walks (i.e., largest connected component in $\mathcal{G}$) is reasonably small, and (d) the minimum required correlation $\varepsilon$ between walks of adjacent sensors is kept low.\vspace{-0mm}

Algorithm~\ref{dadgf} below outlines the key operations of our D$^2$FAS algorithm to be run on each mobile sensor $k$, as detailed previously in Sections~\ref{sec:ddf} and~\ref{sec:dadgf}:\vspace{-3mm}
%The implementation of D$^2$FAS for a team of mobile sensors is shown in . 
%
%\IncMargin{0.5em}
\begin{algorithm}
{\scriptsize
\DontPrintSemicolon
%\SetKwInOut{Input}{inputs}
%  \KwIn{support set $U$, number $K$ of mobile sensors, $k$ (id of current sensor, where $k\in\{1,\cdots,K\}$), $s_k$ (current position of sensor $k$), set $D_k$ of observed road, length $L$ of walk}\BlankLine
\While{true}
{
\tcc{Data fusion (Section~\ref{sec:ddf})}
Construct local summary by (\ref{lsf}) \& (\ref{lsc})\;    
Exchange local summary with every sensor $i\neq k$\; 
Construct global summary by (\ref{gsf}) \& (\ref{gsc})\; 
Predict measurements at unobserved road segments by (\ref{mupitc}) \& (\ref{varpitc})\;
\tcc{Active Sensing (Section~\ref{sec:dadgf})}    
   % Compute the entropy for each $\mathbf{c_k}\in C_k$ and let $\mathbf{c_k^*}$ be the one with maximum entropy\; 
Construct $\Phi_k$ by (\ref{phik})\; 
Exchange $\Phi_k$ with every sensor $i\neq k$\;    
Compute adjacency vector $a_k$ by (\ref{rel2}) \& (\ref{adjv})\;
Exchange adjacency vector with every sensor $i\neq k$\;  
Construct adjacency matrix of coordination graph\;
Find vertex set $\mathcal{V}_n$ of its residing connected component\;
Compute maximum-entropy joint walk $\widehat{w}_{\mathcal{V}_n}$ by (\ref{dmaxent})\;
   % = \arg\max_{w_k} \left| \overline{\Sigma}_{Y_{w_k}Y_{w_k}} \right|$\;
Execute walk $\widehat{w}_k$ and observe its road segments $Y_{\widehat{w}_k}$\;
Update local information $D_k$, $z_{D_k}$, and $s_k$ by (\ref{upDk})\; 
    %, (\ref{upPk}) respectively.
}
\caption{D$^2$FAS$(U, K, L, k, D_k, z_{D_k}, s_k)$}
  \label{dadgf}
}
\end{algorithm}
\vspace{-7mm}
%
\section{Time and Communication Overheads}
\label{overhead}\vspace{-3mm}
In this section, the time and communication overheads of our D$^2$FAS algorithm are analyzed and
compared to that of centralized active sensing (\ref{maxent2}) coupled with the data fusion methods: Full GP (FGP) and SoD (Section~\ref{sec:gpgk}).\vspace{-2mm}
% to be run on each of the $K$ sensors.
%
%
\subsection{Time Complexity}
\label{sec:time}\vspace{-2mm}
%
%Our D$^2$FAS algorithm comprises the data fusion and active sensing components.  
The data fusion component of D$^2$FAS involves computing the local and global summaries and the predictive Gaussian distribution.
%, as shown in Algorithm~\ref{dadgf}.  
To construct the local summary using (\ref{lsf}) and (\ref{lsc}), each sensor has to evaluate $\Sigma_{D_k D_k|U}$ in $\set{O}\hspace{-1mm}\left(|U|^3 + |U| (|D|/K)^2\right)$ time and invert it in $\set{O}\hspace{-1mm}\left((|D|/K)^3\right)$ time, after which the local summary is obtained in $\set{O}\hspace{-1mm}\left(|U|^2|D|/K + |U|(|D|/K)^2\right)$ time. 
The global summary is computed in $\set{O}\hspace{-1mm}\left(|U|^2 K\right)$ by (\ref{gsf}) and (\ref{gsc}). 
Finally, the predictive Gaussian distribution is derived in $\set{O}\hspace{-1mm}\left(|U|^3 + |U||Y|^2 \right)$ time using (\ref{mupitc}) and (\ref{varpitc}).
Supposing $|Y|\leq |U|$ for simplicity, the time complexity of data fusion is then $\set{O}\hspace{-1mm}\left((|D|/K)^3 + |U|^3 + |U|^2 K\right)$.  

%In the second aspect, maximum entropy strategy is computed in decentralized fashion by (\ref{dmaxent}).  
Let the maximum out-degree of $G$ be denoted by $\delta$. Then, each sensor has to consider $\Delta\triangleq\delta^L$ possible walks of length $L$.
The active sensing component of D$^2$FAS involves computing $\Phi_k$ in $\set{O}\hspace{-1mm}\left(\Delta L |U|^2\right)$ time, $a_k$ in $\set{O}\hspace{-1mm}\left(\Delta^{2} L^2 |U| K\right)$ time, its residing connected component in $\set{O}\hspace{-1mm}\left(\kappa^2\right)$ time, and the maximum-entropy joint walk by (\ref{varpitc}) and (\ref{dmaxent}) with the following incurred time: 
The largest connected component of $\kappa$ sensors in $\mathcal{G}$ has to consider $\Delta^{\kappa}$ possible joint walks. Note that 
$\overline{\Sigma}_{Y_{w_{\mathcal{V}_n}} Y_{w_{\mathcal{V}_n}}} = \mbox{diag}\hspace{-1mm}\left((\Sigma_{Y_{w_k}Y_{w_k}|U})_{k\in \mathcal{V}_n}\right) + \Sigma_{Y_{w_{\mathcal{V}_n}}U}\ddot{\Sigma}_{UU}^{-1}\Sigma_{UY_{w_{\mathcal{V}_n}}}$ where $\mbox{diag}(B)$ constructs a diagonal matrix by placing vector $B$ on its diagonal. By exploiting $\Phi_k$, the diagonal and latter matrix terms for all possible joint walks can be computed in $\set{O}\hspace{-1mm}\left(\kappa\Delta(L |U|^2 + L^2 |U|)\right)$ and
$\set{O}\hspace{-1mm}\left(\kappa^2\Delta^{2}L^2 |U|\right)$ time, respectively.
For each joint walk $w_{\mathcal{V}_n}$, evaluating the determinant of $\overline{\Sigma}_{Y_{w_{\mathcal{V}_n}} Y_{w_{\mathcal{V}_n}}}$ incurs $\set{O}\hspace{-1mm}\left((\kappa L)^3 \right)$ time.
%$\set{O}\hspace{-1mm}\left(\kappa L|U|^2 + (\kappa L)^3 \right)$ time.
%the posterior covariance is an $L\times L$ matrix in which each element is evaluated in $\set{O}\hspace{-1mm}\left(|U|^2\right)$ after inversion of .
Therefore, the time complexity of active sensing is $\set{O}\hspace{-1mm}\left(\kappa\Delta L |U|^2 + \Delta^{2} L^2 |U| (K+\kappa^2) +\Delta^{\kappa} (\kappa L)^3\right)$. 

%comparision
Hence, the time complexity of our D$^2$FAS algorithm is $\set{O}\hspace{-0mm}((|D|/K)^3 + |U|^2 (|U| + K + \kappa\Delta L) + \Delta^{2} L^2 |U| (K+\kappa^2)+ \Delta^{\kappa}(\kappa L)^3)$. In contrast, the time incurred by centralized active sensing coupled with FGP and SoD are, respectively, $\set{O}\hspace{-1mm}\left(|D|^3+\Delta^{K}KL(|D|^2 + (KL)^2)\right)$ and $\set{O}\hspace{-1mm}\left(|U|^3|D|+\Delta^{K}KL(|U|^2 + (KL)^2)\right)$. It can be observed that D$^2$FAS can scale better with large $|D|$ (i.e., number of observations) and $K$ (i.e., number of sensors). The scalability of D$^2$FAS vs. FGP and SoD will be further evaluated empirically in Section~\ref{sec:exp}.\vspace{-2mm}
%
%To compare time complexity with FGP and SoD (shown in Table \ref{comp}), considering $K=1$ first, DADGF requires the same $\set{O}\hspace{-1mm}\left(|D|^3\right)$ time as FGP for kernel matrix construction and inversion, but is less efficient than SoD which requires $\set{O}\hspace{-1mm}\left(|U|^3|D|\right)$.  After that, DADGF evaluates one covariance value as efficient as SoD in $\set{O}\hspace{-1mm}\left(|U|^2\right)$ time, while FGP requires $\set{O}\hspace{-1mm}\left(|D|^2\right)$. In this specific case, DADGF is more efficient than FGP and less efficient than SoD.  Nevertheless, when K increases, DADGF can improve the efficiency of the kernel construction due to parallel computation. In contrast, the time for covariance evaluation of SoD boosts exponentially in $K$.  Therefore, in terms of efficiency, DADGF significantly outperforms FGP, and potentially outperforms SoD, especially when the number of observations is limited and the number of mobile sensors is large.
%%%%%%%%%%%%%%%%%%%%%
%
\subsection{Communication Complexity}
\label{sec:comm}\vspace{-2mm}
%
Let the communication overhead be defined as the size of each broadcast message.
%Therefore, the communication complexity is reduced to the size of exchanging message.  
%Recall from Algorithm~\ref{dadgf} (i.e., D$^2$FAS) that, in each iteration, 
Recall from the data fusion component of D$^2$FAS in Algorithm~\ref{dadgf} that, in each iteration, each sensor broadcasts a $\set{O}\hspace{-1mm}\left(|U|^2\right)$-sized summary encapsulating its local observations, which is robust against communication failure. 
In contrast,
%In  (step 2), each mobile sensor exchanges local summary in  size with the other mobile sensors. 
FGP and SoD require each sensor to broadcast, in each iteration, a $\set{O}\hspace{-1mm}\left(|D|/K\right)$-sized message comprising exactly its local observations to handle communication failure.
%have $O(1)$ communication complexity, when exchange only the incoming data points. However, if the communication is not fully reliable which is highly possible in real world, the message may get lost during communication failure. Without exchanging the historical data points, the consistency of the global set $D$ is corrupted. Therefore, the safer way is to exchange all local data which requires $O(|D|)$.  In our approach, the consistency is always enforced by the global summary which is not affected by past communication failure. 
If the number of local observations grows to be larger in size than a local summary of predefined size, then the data fusion component of D$^2$FAS is more scalable than FGP and SoD in terms of communication overhead.
For the partially decentralized active sensing component of D$^2$FAS, each sensor broadcasts $\set{O}\hspace{-1mm}\left(\Delta L|U| \right)$-sized $\Phi_k$ and $\set{O}\hspace{-1mm}\left(K\right)$-sized $a_k$ messages.%\vspace{-3mm}
%
%\begin{table*}
%\centering
%\caption{Comparison of models for active sensing with mobile sensors} 
%\begin{tabular}{|c|c|c|} \hline
%  Model &Time complexity& Communication complexity\\
  % & Coordination Type\\ 
%  \hline
%  FGP   &$\set{O}\hspace{-1mm}\left(|D|^3+\Delta^{KL}(KL|D|^2 + (KL)^3)\right)$ & $\set{O}\hspace{-1mm}\left(|D|\right)$\\ \hline
%  SoD   &$\set{O}\hspace{-1mm}\left(|U|^3|D|+\Delta^{KL}(KL|U|^2 + (KL)^3)\right)$ & $\set{O}\hspace{-1mm}\left(|D|\right)$\\ \hline
%  D$^2$FAS & $\set{O}\hspace{-1mm}\left((|D|/K)^3+|U|^3+|U|^2 K + \Delta^L(L|U|^2 +L^3)\right)$             &$\set{O}\hspace{-1mm}\left(|U|^2\right)$\\ \hline
%\end{tabular}
%\label{comp}
%\end{table*}
%
%pro & con
\section{Experiments and Discussion}
\label{sec:exp}\vspace{-3mm}
%
\begin{wrapfigure}{l}{4.1cm}
%\centering
%\setlength\fboxrule{0.1pt}
\vspace{-4mm}\includegraphics[scale=0.40]{dom775.pdf}\vspace{-4mm}
%\caption{}
%\caption{Traffic phenomenon (i.e., speeds (km/h) of road segments) over an urban road network in Anonymous area, Anonymity during evening peak hours on April $20$, $2011$. It comprises $775$ road segments including highways, arterials, slip roads, etc. The mean speed is $48.8$ km/h and the population standard deviation is $20.5$ km/h.}
%\label{dom775}\vspace{-4mm}
\end{wrapfigure}
This section evaluates the predictive performance, time efficiency, and scalability of our D$^2$FAS algorithm on a real-world traffic phenomenon (i.e., speeds (km/h) of road segments) over an urban road network
%, as shown and detailed in Fig.~\ref{dom775}.
(top figure) in Tampines area, Singapore during evening peak hours on April $20$, $2011$.
It comprises $775$ road segments including highways, arterials, slip roads, etc. The mean speed is $48.8$ km/h and the standard deviation is $20.5$ km/h.
%The road network consists of $775$ road segments including both highway and arterial segments.

The performance of D$^2$FAS is compared to that of centralized active sensing (\ref{maxent2}) coupled with the state-of-the-art data fusion methods: full GP (FGP) and SoD (Section~\ref{sec:gpgk}).  
%Additionally, the velocity on each segment is averaged over measurements logged by passing bus at 5-minute interval during 5 P.M. to 7 P.M., on April 20, 2011. 
%assumptions
A network of $K$ mobile sensors is tasked to explore the road network to gather a total of up to $960$ observations. 
To reduce computational time, each sensor repeatedly computes and executes maximum-entropy walks of length $L=2$ (instead of computing a very long walk), unless otherwise stated. For D$^2$FAS and SoD, $|U|$ is set to $64$ .  
For the active sensing component of D$^2$FAS, $\varepsilon$ is set to $0.1$, unless otherwise stated.
The experiments are run on a Linux PC with Intel$\circledR$ Core$^{\mbox{\tiny{TM}}}2$ Quad CPU Q$9550$ at $2.83$~GHz.\vspace{-3mm}
%\begin{figure}
%\centering
%\setlength\fboxrule{0.1pt}
%\fbox{\includegraphics[scale=0.68]{dom775}}\vspace{-3mm}
%\caption{Traffic phenomenon (i.e., speeds (km/h) of road segments) over an urban road network in Anonymous area, Anonymity during evening peak hours on April $20$, $2011$. It comprises $775$ road segments including highways, arterials, slip roads, etc. The mean speed is $48.8$ km/h and the population standard deviation is $20.5$ km/h.}
%\label{dom775}\vspace{-4mm}
%\end{figure}
%
%
\subsection{Performance Metrics}\vspace{-3mm}
%
%The performance of the algorithms is evaluated in terms of prediction accuracy, run-time efficiency and scalability. 
The first metric evaluates the predictive performance of a tested algorithm: It measures the \emph{root mean squared error} (RMSE)
%
%\begin{equation*}
$\sqrt{|V|^{-1}\sum_{s\in V}\left(z_s-\widehat{\mu}_s\right)^2}$
%\label{eq:rmse}
%\end{equation*}
%
over the entire domain $V$ of the road network that is incurred by the predictive mean $\widehat{\mu}_s$ of the tested algorithm, specifically, using (\ref{pmu}) of FGP,
(\ref{eq:fsod}) of SoD, or (\ref{mupitc}) of D$^2$FAS.
%resulting from using the speed measurements at the set $D$ of observed road segments and the 
%where $\mu_{s|D}$ is the predictive velocity value. 
The second metric evaluates the time efficiency and scalability of a tested algorithm by measuring its incurred time; for D$^2$FAS, the maximum of the time incurred by all subsets $\mathcal{V}_1,\ldots,\mathcal{V}_{\mathcal{K}}$ of sensors is recorded.\vspace{-3mm}
%The run-time efficiency and scalability are evaluated by comparing the incurred time from the experiment.  
%
%
\subsection{Results and Analysis}
\label{ppte}\vspace{-3mm}
\begin{figure}[h]
\begin{tabular}{ccc}
\hspace{-2mm}\includegraphics[scale=0.143]{acc4.pdf} & \hspace{-4.5mm}\includegraphics[scale=0.143]{acc6.pdf} &  \hspace{-4.5mm}\includegraphics[scale=0.143]{acc8.pdf}\vspace{-1mm}\\
\hspace{-2mm}(a) $K = 4$ & \hspace{-4.5mm}(b) $K = 6$ &  \hspace{-4.5mm}(c) $K = 8$\\ 
 \hspace{-2mm}\includegraphics[scale=0.143]{eff4.pdf} & \hspace{-4.5mm}\includegraphics[scale=0.143]{eff6.pdf} &  \hspace{-4.5mm}\includegraphics[scale=0.143]{eff8.pdf}\vspace{-1mm}\\
\hspace{-2mm}(d) $K = 4$ & \hspace{-4.5mm}(e) $K = 6$ &  \hspace{-4.5mm}(f) $K = 8$\vspace{-3mm} 
\end{tabular}
\caption{Graphs of (a-c) predictive performance and (d-f) time efficiency vs. total no. $|D|$ of observations gathered by varying number $K$ of mobile sensors.}
\label{fig3}\vspace{-2mm}
\end{figure}
{\bf Predictive performance and time efficiency.} 
Fig.~\ref{fig3} shows results of the performance of the tested algorithms averaged over $40$ randomly generated starting sensor locations with varying number $K=4, 6, 8$ of sensors. It can be observed that D$^2$FAS is significantly more time-efficient 
and scales better with increasing number $|D|$ of observations 
(Figs.~\ref{fig3}d to~\ref{fig3}f) while achieving predictive performance close to that of centralized active sensing coupled with FGP and SoD (Figs.~\ref{fig3}a to~\ref{fig3}c).
Specifically, D$^2$FAS is about $1, 2, 4$ orders of magnitude faster than centralized active sensing coupled with FGP and SoD for $K = 4,6,8$ sensors, respectively.
%This indicates that the capability of D$^2$FAS is sophisticated enough for modeling and predicting the velocity phenomena on road network. 
%Additionally, it implies the conditional independence assumption for decentralizing joint entropy computing is reasonable.
%Fig.~\ref{eff4} compares the performance of three approaches, with 4 mobile sensors exploring the road network.  
%The team of mobile sensors will gather a total number of 960 observations in the velocity field over road network (240 observations/sensor) with planning horizon 2. 
%finding -> explanation
%In terms of incurred time, Fig. \ref{eff4}a shows that D$^2$FAS significantly outperforms FGP and SoD. This result is expected from the time complexity comparison of three methods in Section \ref{sec:time}.
%We find that D$^2$FAS can outperform SoD when number of sensors is four
%with less than 960 observations. We believe that, in practice, the team
%size of mobile sensors is highly likely to exceed four, thus D$^2$FAS
%is more favorable than SoD. 
%Fig. \ref{eff4}b shows that D$^2$FAS has nearly the same predictive error as FGP and SoD. 
%eff4
%With the same experimental setting except for using 6 mobile sensors (160 observations/sensor), Fig. \ref{eff6} shows very similar results. In addition, we can observe that a sensor with FGP requires around 120 seconds to plan a path of length 2, which is impractical for online application.  
%sum-up
%To sum up, our D$^2$FAS algorithm significantly outperform the FGP and SoD approach, and achieve comparable prediction accuracy. 
%
%\subsection{Scalability}
%\label{sec:scale}
%
\begin{figure}
\begin{tabular}{ccc}
\hspace{-2.5mm}\includegraphics[scale=0.143]{scale_d2fas_err.pdf} & \hspace{-4mm}\includegraphics[scale=0.143]{scale_fgp_err.pdf} &  \hspace{-4mm}\includegraphics[scale=0.143]{scale_sod_err.pdf}\vspace{-1mm}\\
\hspace{-2.5mm}(a) D$^2$FAS & \hspace{-4mm}(b) FGP &  \hspace{-4mm}(c) SoD\\ 
 \hspace{-2.5mm}\includegraphics[scale=0.143]{scale_d2fas_tm.pdf} & \hspace{-4mm}\includegraphics[scale=0.143]{scale_fgp_tm.pdf} &  \hspace{-4mm}\includegraphics[scale=0.143]{scale_sod_tm.pdf}\vspace{-1mm}\\
\hspace{-2.5mm}(d) D$^2$FAS & \hspace{-4mm}(e) FGP &  \hspace{-4mm}(f) SoD\vspace{-3mm} 
\end{tabular}
\caption{Graphs of (a-c) predictive performance and (d-f) time efficiency vs. total no. $|D|$ of observations gathered by varying number $K$ of mobile sensors.}
\label{fig6}\vspace{-4mm}
\end{figure}

{\bf Scalability of D$^2$FAS.} Using the same results as that in Fig.~\ref{fig3}, Fig.~\ref{fig6} plots them differently to reveal the scalability of the tested algorithms with increasing number $K$ of sensors.
Additionally, we provide results of the performance of D$^2$FAS for $K=10, 20, 30$ sensors; 
such results are not available for centralized active sensing coupled with FGP and SoD due to extremely long incurred time.
It can be observed from Figs.~\ref{fig6}a to~\ref{fig6}c that the predictive performance of all tested algorithms improve with a larger number of sensors because each sensor needs to execute fewer walks and its performance is therefore less adversely affected by its myopic selection (i.e., $L=2$) of maximum-entropy walks. As a result, more informative unobserved road segments are explored.
%From Figs. \ref{scale1}d-f , we observe the same trend that the prediction error is reduced as K increases. This is because, with a larger team of mobile sensors, the more informative segments could be explored with higher chance. 

As shown in Fig.~\ref{fig6}d, when the randomly placed sensors gather their initial observations (i.e., $|D|<400$), the time incurred by D$^2$FAS is higher for greater $K$ due to larger subsets of sensors being formed to coordinate their walks (i.e., larger $\kappa$). As more observations are gathered (i.e., $|D|\geq 400$), its partially decentralized active sensing component directs the sensors to explore further apart from each other in order to maximize the entropy of their walks.
This consequently decreases $\kappa$, leading to a reduction in incurred time.
Furthermore, as $K$ increases from $4$ to $20$, the incurred time decreases due to its decentralized data fusion component that can distribute the computational load among a greater number of sensors.
When the road network becomes more crowded from $K=20$ to $K=30$ sensors, the incurred time increases slightly due to slightly larger $\kappa$.
%As shown in Fig.~\ref{fig6}d, the time incurred by D$^2$FAS decreases due to its decentralized data fusion component that can distribute the computational load among a greater number of sensors.
In contrast, Figs.~\ref{fig6}e and~\ref{fig6}f show that the time taken by FGP and SoD increases significantly primarily due to their centralized active sensing incurring exponential time in $K$.
Hence, the scalability of our D$^2$FAS algorithm in the number of sensors allows the deployment of a larger-scale mobile sensor network (i.e., $K\geq 10$) to achieve more accurate traffic modeling and prediction (Figs.~\ref{fig6}a to~\ref{fig6}c).
\begin{figure}
\begin{tabular}{ccc}
 \hspace{-2mm}\includegraphics[scale=0.143]{scale_pitc_tm.pdf} & \hspace{-4mm}\includegraphics[scale=0.143]{dscal_fgp_tm.pdf} &  \hspace{-4mm}\includegraphics[scale=0.143]{dscal_sod_tm.pdf}\vspace{-1mm}\\
\hspace{-2mm}(a) D$^2$FAS & \hspace{-4mm}(b) FGP &  \hspace{-4mm}(c) SoD\vspace{-3mm} 
\end{tabular}
\caption{Graphs of time efficiency vs. total no. $|D|$ of observations gathered by varying number $K$ of sensors.}
\label{fig6b}\vspace{-5mm}
\end{figure}

{\bf Scalability of data fusion.} Fig.~\ref{fig6b} shows results of the scalability of the tested data fusion methods with increasing number $K$ of sensors. In order to produce meaningful results for fair comparison, the same active sensing component has to be coupled with the data fusion methods and its incurred time kept to a minimum.
As such, we impose the use of a fully decentralized active sensing component to be performed by each mobile sensor $k$:
%\vspace{-1mm} 
%\begin{equation}
%\begin{array}{c}
%\hspace{-3mm}
$w^{\ast}_k = \mathop{\arg\max}_{w_k} \mathbb{H}[Z_{Y_{w_k}}| Z_{D}]$.
%\vspace{-1mm}
%\\
 % \displaystyle\widehat{\mathbb{H}}[Z_{Y_{w_k}}]\triangleq \frac{1}{2}\log(2\pi e)^{\left|Y_{w_k}\right|}\left|\widehat{\Sigma}_{Y_{w_k}Y_{w_k}}\right|
 % \end{array}
%\label{dmaxent3}
%\end{equation}
For D$^2$FAS, this corresponds exactly to (\ref{dmaxent}) by setting a large enough $\varepsilon$ (in our experiments, $\varepsilon=2$) to yield $\kappa=1$; consequently, 
%its active sensing component is simplified to computing the maximum-entropy walk; 
computational and communicational operations pertaining to the coordination graph can be omitted.

It can be seen from Fig.~\ref{fig6b}a that the time incurred by the decentralized data fusion component of D$^2$FAS decreases with increasing $K$, as explained previously.
In contrast, the time incurred by FGP and SoD increases (Fig.~\ref{fig6b}b and~\ref{fig6b}c):
As discussed above, a larger number of sensors results in a greater quantity of more informative unique observations to be gathered (i.e., fewer repeated observations), which increases the time needed for data fusion. 
When $K\geq 10$, D$^2$FAS is at least $1$ order of magnitude faster than FGP and SoD.
It can also be observed that D$^2$FAS scales better with increasing number of observations.
So, the real-time performance and scalability 
of D$^2$FAS's decentralized data fusion enable it to be used for persistent large-scale traffic modeling and prediction where a large number of observations and sensors (including static and passive ones) are expected to be available.
%Fig. \ref{scale1} shows how the run-time and prediction error are affected by scaling up the amount of mobile sensors.  We simply modify the previous experimental setting by varying the number of mobile sensors (K=2,4,6,8,10,20,30).  From Fig. \ref{scale1}a, we find that the incurred time decreases when mobile sensors scale up, while SoD and FGP show the opposite trend in Figs. \ref{scale1}b-c respectively. The opposite trend between decentralized method and centralized methods is well explained by Table \ref{comp}. Clearly, when K scales up, the time complexity of D$^2$FAS is reduced due to parallel computation of the decentralized data fusion, while FGP and SoD have to evaluate the joint entropy with exponentially time in K.  Due to the same reason, we can observe that D$^2$FAS can scale up to 30 sensors in extremely short time, in contrast, SoD can scale up to 8 sensors and FGP, even worse, can only scale up to 6 sensors.
%
%The approach with better scalability is more favorable, as higher accuracy can achieved by deploying more mobile sensors. This argument is supported by the results that, as we can see, the lowest prediction error is around 6 with 30 sensors (shown in Fig.~\ref{scale1}d), however the one with less sensors is larger that 8 (shown in Figs \ref{scale1}e,\ref{scale1}f).
%\begin{figure}
%\begin{tabular}{cc} 
%\hspace{-2mm}\includegraphics[scale=0.22]{fig6/dscal_pitc_err.eps} &
%\hspace{-3mm}\includegraphics[scale=0.22]{fig6/dscal_pitc_tm.eps} \\ 
% \hspace{-2mm} (a) D$^2$FAS & \hspace{-3mm}(d) D$^2$FAS\\
%\hspace{-2mm}\includegraphics[scale=0.22]{fig6/dscal_fgp_err.eps} &
%\hspace{-3mm}\includegraphics[scale=0.22]{fig6/dscal_fgp_tm.eps} \\ 
%\hspace{-2mm}  (b) FGP & \hspace{-3mm}(e) FGP\\
%\hspace{-2mm}\includegraphics[scale=0.22]{fig6/dscal_sod_err.eps} &
%\hspace{-3mm}\includegraphics[scale=0.22]{fig6/dscal_sod_tm.eps} \\ 
% \hspace{-2mm} (c) SoD & \hspace{-3mm}(f) SoD
%\end{tabular}
%\caption{Graphs of (a-c) predictive performance and (d-f) time efficiency vs. total no. $|D|$ of observations gathered by varying number $K$ of sensors.}
%\label{fig6}
%\end{figure}
%
%\subsection{Planning Horizon}
\begin{figure}
\begin{tabular}{ccc}
\hspace{-2.5mm}\includegraphics[scale=0.143]{hor_d2fas_err.pdf} 
%\hspace{-2.5mm}\includegraphics[scale=0.143]{fig7/hor_pitc_err.eps} 
& \hspace{-4mm}\includegraphics[scale=0.143]{hor_fgp_err.pdf} &  \hspace{-4mm}\includegraphics[scale=0.143]{hor_sod_err.pdf}\vspace{-1mm}\\
\hspace{-2.5mm}(a) D$^2$FAS & \hspace{-4mm}(b) FGP &  \hspace{-4mm}(c) SoD\\ 
%\hspace{-2.5mm}\includegraphics[scale=0.143]{fig7/hor_pitc_err.eps} 
\hspace{-2.5mm}\includegraphics[scale=0.143]{hor_d2fas_tm.pdf} 
% \hspace{-2.5mm}\includegraphics[scale=0.143]{fig7/hor_pitc_tm.eps} 
& \hspace{-4mm}\includegraphics[scale=0.143]{hor_fgp_tm.pdf} &  \hspace{-4mm}\includegraphics[scale=0.143]{hor_sod_tm.pdf}\vspace{-1mm}\\
\hspace{-2.5mm}(d) D$^2$FAS & \hspace{-4mm}(e) FGP &  \hspace{-4mm}(f) SoD\vspace{-3mm} 
\end{tabular}
\caption{Graphs of (a-c) predictive performance and (d-f) time efficiency vs. total no. $|D|$ of observations gathered by $2$ mobile sensors with varying length $L$ of maximum-entropy joint walks.}
\label{fig5}\vspace{-3mm}
\end{figure}

{\bf Varying length $L$ of walk.} Fig.~\ref{fig5} shows results of the performance of the tested algorithms with varying length $L=2, 4, 6, 8$ of maximum-entropy joint walks; we choose to experiment with just $2$ sensors since Figs.~\ref{fig6} and~\ref{fig6b} reveal that a smaller number of sensors produce poorer predictive performance and higher incurred time with large number of observations for D$^2$FAS. It can be observed that the predictive performance of all tested algorithms improve with increasing walk length $L$ because the selection of maximum-entropy joint walks is less myopic.
The time incurred by D$^2$FAS increases due to larger $\kappa$ but grows more slowly and is lower than that incurred by centralized active sensing coupled with FGP and SoD.
Specifically, when $L=8$, D$^2$FAS is at least $1$ order of magnitude faster (i.e., average of $60$~s) than centralized active sensing coupled with SoD (i.e., average of $>732$~s) and FGP (i.e., not available due to excessive incurred time).
%In particular, D$^2$FAS can select walks of length $L=8$ to yield better  
%When $L$ increases to $10$, the incurred time increases to about $10$ seconds, which is reasonable in practice. 
%L=2   FGP: 21=0.0875   SoD: 48=0.2   D2FAS: 14=0.058333
%L=4   FGP: 342=2.85   SoD: 44=0.366667   D2FAS: 17=0.141667
%L=6   FGP: 29919s=373.9875   SoD:840s=10.5   D2FAS:221s=2.7625
%L=8   FGP: na   SoD:43926s=732.1   D2FAS:3580s=59.666667
% I think maybe we should divide by the number of iterations.
% For example, to collect 960 observations, L=2 : 240 iterations (K=2), L= 4 : 120 iterations
% L=6 : 80 iterations, L=8 : 60 iterations
Also, notice from Figs.~\ref{fig6}a and~\ref{fig6}d that if a large number of sensors (i.e., $K=30$) is available, D$^2$FAS can select shorter walks of $L=2$ to be significantly more time-efficient (i.e., average of $>3$ orders of magnitude faster) while achieving predictive performance comparable to that of SoD with $L=8$ and FGP with $L=6$.\vspace{-2mm}

%By deploying a larger number of sensors, the incurred time is expected to decrease while improving the predictive performance.\vspace{-2mm}
%\begin{figure}
%\centering
%\begin{tabular}{cc} 
%\hspace{-2mm}\includegraphics[scale=0.22]{fig5/hor_pitc_err.eps} &
%\hspace{-3mm}\includegraphics[scale=0.22]{fig5/hor_pitc_tm.eps} \\ 
%\hspace{-2mm}(a) & \hspace{-3mm}(b)
%\end{tabular}
%\caption{Graphs of (a) predictive performance and (b) time efficiency vs. total no. $|D|$ of observations gathered by $2$ mobile sensors running D$^2$FAS with varying length $L$ of maximum-entropy walks.}
%\label{fig5}
%\end{figure}
%
\section{Conclusion}\vspace{-3mm}
%
This paper describes a decentralized data fusion and active sensing algorithm for modeling and predicting spatiotemporal traffic phenomena with mobile sensors.
Analytical and empirical results have shown that our D$^2$FAS algorithm is significantly more time-efficient and scales better with increasing number of observations and sensors while achieving predictive performance close to that of state-of-the-art centralized active sensing coupled with FGP and SoD.
Hence, D$^2$FAS is practical for deployment in a large-scale mobile sensor network to achieve persistent and accurate traffic modeling and prediction.
For our future work, we will assume that each sensor can only communicate locally with its neighbors (instead of assuming all-to-all communication) and develop a \emph{distributed} data fusion approach to efficient and scalable approximate GP prediction based on D$^2$FAS and consensus filters \citep{saber05a}.
%The DADGF algorithm is remarkably efficient, since the time complexity of computing joint entropy is independent ofthe number of mobile sensors and local data are assimilated faster as the number of mobile sensors increase. The team of mobile sensors employing DADGF algorithm is highly scalable as communication complexity is independent of the increasing data. Additionally, the capability of DADGF is proved to be equivalent to a sophisticated centralized sparse GP model (PITC).  Empirical evaluation on real-world urban road network shows, with a support set of reasonable size, that DADGF significantly outperforms the state-of-art centralized GP-based active sensing models in terms of run-time efficiency and scalability, and achieves comparable prediction accuracy.

{\bf Acknowledgments.} This work was supported by Singapore-MIT Alliance Research and Technology (SMART) Subaward Agreement $14$ R-$252$-$000$-$466$-$592$.

\clearpage

\bibliographystyle{natbib}
\bibliography{mybib}
%% Use plainnat to work nicely with natbib. 

\ifthenelse{\value{sol}=1}{
\clearpage

\appendix
%\section{Matrix Inversion Lemma} 
%Let $Z$ be a $n\times n$ matrix, $W$ be a $m\times m$ matrix, and $U$ and $V$ be $n\times m$ matrices. Supposing the matrix inverses of $Z$ and $W$ exist, 
%
%\begin{equation} 
 % \begin{array}{l}
%(Z+UWV^T)^{-1}\\ 
%= Z^{-1}-Z^{-1}U(W^{-1}+V^TZ^{-1}U)^{-1}V^TZ^{-1} \ .
%\end{array} 
 % \label{eq:matinv} 
%\end{equation} 
% 
%
\section{Proof of Theorem~\ref{thmdad}B}
\label{equ}
%We will prove that the PITC model is equivalent to the  $D^2GF$ model.   
%The centralized partially independent training conditional (PITC) approximation \citep{candela05} yields the predictive Gaussian distribution $\mathcal{N}(\mu_Y^{\mbox{\scriptsize PITC}},\Sigma_{YY}^{\mbox{\scriptsize PITC}})$ of the measurements at any set $Y$ of unobserved road segments where
%\begin{equation} 
%  \mu_Y^{\mbox{\scriptsize PITC}}\triangleq\mu_Y+\Gamma_{YD}\left(\Gamma_{DD}+\Lambda\right)^{-1}\left(z_D-\mu_D \right)
%  \label{eq:pitc.mu} 
%\end{equation}
%\begin{equation}
 % \Sigma_{YY}^{\mbox{\scriptsize PITC}}\triangleq\Sigma_{YY}-\Gamma_{YD}\left(\Gamma_{DD}+\Lambda\right)^{-1}\Gamma_{DY}
  %\label{eq:pitc.var} 
%\end{equation} 
%such that 
%\begin{equation} 
 % Q_{AB} \triangleq
  %\Sigma_{AU}\Sigma_{UU}^{-1}\Sigma_{UB} 
  %\label{eq:sparseq}
%\end{equation} 
%and 
%$U$ is the support set for sparse GP models.
%$\Lambda\triangleq
%\mbox{blockdiag}[\Sigma_{DD|U}]=
%\mbox{blockdiag}[\Sigma_{DD}-\Gamma_{DD}]$ is a $|D|\times |D|$ block diagonal matrix comprising $K$ blocks, each of which is a $|D_k|\times |D_k|$ matrix $\Lambda_k \triangleq\Sigma_{D_k D_k|U}=\Sigma_{D_k D_k}-Q_{D_k D_k}$ for $k=1,\ldots,K$.
%
We have to first simplify the $\Gamma_{YD}\left(\Gamma_{DD}+\Lambda\right)^{-1}$ term in the expressions of $\mu_{Y|D}^{\mbox{\tiny PITC}}$ (\ref{eq:pitc.mu}) and $\Sigma_{YY|D}^{\mbox{\tiny PITC}} $ (\ref{eq:pitc.var}).
  \begin{equation} 
  \hspace{-2mm}
    \begin{array}{l} 
\left(\Gamma_{DD}+\Lambda\right)^{-1}\\
\displaystyle = \left(\Sigma_{DU}\Sigma_{UU}^{-1}\Sigma_{UD}+\Lambda\right)^{-1}\\
\displaystyle = \Lambda^{-1}-\Lambda^{-1}\Sigma_{DU}
    \left(\Sigma_{UU}+\Sigma_{UD}\Lambda^{-1}\Sigma_{DU}\right)^{-1}
    \Sigma_{UD}\Lambda^{-1}\\ 
%\displaystyle = \Lambda^{-1}-\Lambda^{-1}\Sigma_{DU}\hspace{-1mm}\left(\hspace{-1mm}\Sigma_{UU}+\sum_{k=1}^{K}\Sigma_{UD_k}\Lambda_k^{-1}\Sigma_{D_kU}\hspace{-1mm}\right)^{\hspace{-2mm}-1}\hspace{-1mm}\Sigma_{UD}\Lambda^{-1}\\ 
%\displaystyle = \Lambda^{-1}-\Lambda^{-1}\Sigma_{DU}\left(\Sigma_{UU}+\sum_{k=1}^{K}\dot{\Sigma}_{UU}^k \right)^{-1}\Sigma_{UD}\Lambda^{-1}\\
\displaystyle =\Lambda^{-1}-\Lambda^{-1}\Sigma_{DU}\ddot{\Sigma}_{UU}^{-1}\Sigma_{UD}\Lambda^{-1} \ .
  \end{array} 
    \label{eq:pitck} 
\end{equation}
The second equality follows from matrix inversion lemma.
The last equality is due to 
\begin{equation}
    \begin{array}{l}
\displaystyle\Sigma_{UU}+\Sigma_{UD}\Lambda^{-1}\Sigma_{DU}\\
=\displaystyle \Sigma_{UU}+\sum_{k=1}^{K}\Sigma_{UD_k}\Sigma_{D_k D_k|U}^{-1}\Sigma_{D_kU}\\
= \displaystyle \Sigma_{UU}+\sum_{k=1}^{K}\dot{\Sigma}_{UU}^k =\displaystyle \ddot{\Sigma}_{UU} \ .
  \end{array} 
    \label{eq:aa} 
\end{equation}
Using (\ref{eq:sparseq}) and (\ref{eq:pitck}),
  \begin{equation} 
  %\hspace{-4mm}
    \begin{array}{l} 
\Gamma_{YD}\left(\Gamma_{DD}+\Lambda\right)^{-1}\\
=\displaystyle\Sigma_{YU}\Sigma_{UU}^{-1}\Sigma_{UD}\left(\Lambda^{-1}-\Lambda^{-1}\Sigma_{DU}\ddot{\Sigma}_{UU}^{-1}\Sigma_{UD}\Lambda^{-1}\right)\\
=\displaystyle\Sigma_{YU}\Sigma_{UU}^{-1}\left(\ddot{\Sigma}_{UU}-\Sigma_{UD}\Lambda^{-1}\Sigma_{DU}\right)\ddot{\Sigma}_{UU}^{-1}\Sigma_{UD}\Lambda^{-1}\\
=\displaystyle\Sigma_{YU}\ddot{\Sigma}_{UU}^{-1}\Sigma_{UD}\Lambda^{-1}
  \end{array} 
    \label{eq:pitck2} 
\end{equation}
The third equality is due to (\ref{eq:aa}).

From (\ref{eq:pitc.mu}),
%Then, we can rewrite (\ref{eq:pitc.mu}) as 
%
\begin{align*} 
%  \hspace{-2mm}
 \mu_{Y|D}^{\mbox{\tiny PITC}}
&=\displaystyle\mu_Y+\Gamma_{YD}\left(\Gamma_{DD}+\Lambda\right)^{-1}\left(z_D-\mu_D \right)\\
%\displaystyle = \mu_Y+\Sigma_{YU}\Sigma_{UU}^{-1}\Sigma_{UD}\left(\Lambda^{-1}-\Lambda^{-1}\Sigma_{DU}\ddot{\Sigma}_{UU}^{-1}\Sigma_{UD}\Lambda^{-1}\right)\left(z_D-\mu_D \right)\\
%\displaystyle = \mu_Y+\Sigma_{YU}\Sigma_{UU}^{-1}\left(\ddot{\Sigma}_{UU}-\Sigma_{UD}\Lambda^{-1}\Sigma_{DU}\right)\ddot{\Sigma}_{UU}^{-1}\Sigma_{UD}\Lambda^{-1}\left(z_D-\mu_D \right)\\
&=\displaystyle \mu_Y+\Sigma_{YU}\ddot{\Sigma}_{UU}^{-1}\Sigma_{UD}\Lambda^{-1}\left(z_D-\mu_D \right)\\
&=\displaystyle \mu_Y+\Sigma_{YU}\ddot{\Sigma}_{UU}^{-1}\ddot{z}_U\\
&= \overline{\mu}_{Y} \ .
  \label{eq:pitc.mur}
\end{align*}
The second equality is due to (\ref{eq:pitck2}).
The third equality follows from 
$
\Sigma_{UD}\Lambda^{-1}\left(z_D-\mu_D \right)
= \sum_{k=1}^{K}\Sigma_{UD_k}\Sigma_{D_k D_k|U}^{-1}\left(z_{D_k}-\mu_{D_k} \right) = \sum_{k=1}^{K} \dot{z}^k_U =\ddot{z}_U$.

From (\ref{eq:pitc.var}),
%we can rewrite (\ref{eq:pitc.var}) as 
%
\begin{equation*}
  \begin{array}{l} 
  \Sigma_{YY|D}^{\mbox{\tiny PITC}}\\
   \displaystyle=\Sigma_{YY}-\Gamma_{YD}\left(\Gamma_{DD}+\Lambda\right)^{-1}\Gamma_{DY}\\
 \displaystyle    =\Sigma_{YY}-\Sigma_{YU}\ddot{\Sigma}_{UU}^{-1}\Sigma_{UD}\Lambda^{-1}\Sigma_{DU}\Sigma_{UU}^{-1}\Sigma_{UY}\\
 \displaystyle    =\Sigma_{YY}-\Big(\Sigma_{YU}\ddot{\Sigma}_{UU}^{-1}\Sigma_{UD}\Lambda^{-1}\Sigma_{DU}\Sigma_{UU}^{-1}\Sigma_{UY}\\
 \displaystyle \ \ \ \  -\Sigma_{YU}\Sigma_{UU}^{-1}\Sigma_{UY}\Big)
    -\Sigma_{YU}\Sigma_{UU}^{-1}\Sigma_{UY}\\
  \displaystyle   =\Sigma_{YY}- \Sigma_{YU}\ddot{\Sigma}_{UU}^{-1}\left(\Sigma_{UD}\Lambda^{-1}\Sigma_{DU}
    -\ddot{\Sigma}_{UU}\right)\Sigma_{UU}^{-1}\Sigma_{UY}\\
  \displaystyle \ \ \ \   -\Sigma_{YU}\Sigma_{UU}^{-1}\Sigma_{UY}\\
    \displaystyle =\Sigma_{YY}-\left(\Sigma_{YU}\Sigma_{UU}^{-1}\Sigma_{UY}-\Sigma_{YU}\ddot{\Sigma}_{UU}^{-1}\Sigma_{UY}\right)\\
    \displaystyle =\Sigma_{YY}-\Sigma_{YU}\left(\Sigma_{UU}^{-1}-\ddot{\Sigma}_{UU}^{-1}\right)\Sigma_{UY}\\
    = \overline{\Sigma}_{YY} \ .
  \end{array} 
      \label{eq:pitc.varr} 
\end{equation*}
The second equality follows from (\ref{eq:sparseq}) and (\ref{eq:pitck2}).
The fifth equality is due to (\ref{eq:aa}).
 %
\section{Proof of Theorem~\ref{thm2}} 
\label{equ2}
Let $\widetilde{\Sigma}_{Y_{w} Y_{w}} \triangleq \overline{\Sigma}_{Y_{w} Y_{w}} - \widehat{\Sigma}_{Y_{w} Y_{w}}$
%the eigenvalues of $\left(\widehat{\Sigma}_{Y_{w} Y_{w}}\right)^{\hspace{-1mm}-1} \widetilde{\Sigma}_{Y_{w} Y_{w}}$ be denoted by $\lambda_1,\ldots,\lambda_{|Y_{w}|}$,
and $\rho_w$ be the spectral radius of $\left(\widehat{\Sigma}_{Y_{w} Y_{w}}\right)^{\hspace{-1mm}-1} \widetilde{\Sigma}_{Y_{w} Y_{w}}$.
%(i.e., $\rho = \max\hspace{-0.5mm}\left( \left|\lambda_1\right|,\ldots,\left|\lambda_{|Y_{w}|}\right|\right)$.
We have to first bound $\rho_w$ from above.

For any joint walk $w$, $\left(\widehat{\Sigma}_{Y_{w} Y_{w}}\right)^{\hspace{-1mm}-1} \widetilde{\Sigma}_{Y_{w} Y_{w}}$ comprises diagonal blocks of size $\left|Y_{w_{\mathcal{V}_n}}\right|\times \left|Y_{w_{\mathcal{V}_n}}\right|$ with components of value $0$ for $n=1,\ldots,\mathcal{K}$ and off-diagonal blocks of the form $\left(\overline{\Sigma}_{Y_{w_{\mathcal{V}_n}} Y_{w_{\mathcal{V}_n}}}\right)^{\hspace{-1mm}-1}\hspace{-1mm}\overline{\Sigma}_{Y_{w_{\mathcal{V}_n}} Y_{w_{\mathcal{V}_{n'}}}}$ for $n,n' =1,\ldots,\mathcal{K}$ and $n\neq n'$.
We know that any pair of sensors $k\in \mathcal{V}_n$ and $k'\in \mathcal{V}_{n'}$ reside in different connected components of coordination graph $\mathcal{G}$ and are therefore not adjacent.
So, by Definition~\ref{adj},
\begin{equation}
\displaystyle\max_{i,i'}\left|\left[\overline{\Sigma}_{Y_{w_{\mathcal{V}_n}} Y_{w_{\mathcal{V}_{n'}}}}\right]_{ii'}\right|\leq\varepsilon
\label{max2}
\end{equation}
for $n,n' =1,\ldots,\mathcal{K}$ and $n\neq n'$.
Using (\ref{max1}) and (\ref{max2}), each component in any off-diagonal block of $\left(\widehat{\Sigma}_{Y_{w} Y_{w}}\right)^{\hspace{-1mm}-1} \widetilde{\Sigma}_{Y_{w} Y_{w}}$ can be bounded as follows:
\begin{equation}
\displaystyle\max_{i,i'}\left|\left[\left(\overline{\Sigma}_{Y_{w_{\mathcal{V}_n}} Y_{w_{\mathcal{V}_n}}}\right)^{\hspace{-1mm}-1}\hspace{-1mm}\overline{\Sigma}_{Y_{w_{\mathcal{V}_n}} Y_{w_{\mathcal{V}_{n'}}}} \right]_{ii'}\right|\leq\left|Y_{w_{\mathcal{V}_n}}\right|\xi\varepsilon
\label{max3}
\end{equation}
for $n,n' =1,\ldots,\mathcal{K}$ and $n\neq n'$.
It follows from (\ref{max3}) that 
\begin{equation}
\displaystyle\max_{i,i'}\left|\left[\left(\widehat{\Sigma}_{Y_{w} Y_{w}}\right)^{\hspace{-1mm}-1} \widetilde{\Sigma}_{Y_{w} Y_{w}}\right]_{ii'}\right|\leq\max_{n}\left|Y_{w_{\mathcal{V}_n}}\right|\xi\varepsilon\leq L\kappa\xi\varepsilon\ .
\label{max4}
\end{equation}
The last inequality is due to
$\displaystyle\max_{n}\left|Y_{w_{\mathcal{V}_n}}\right|\leq L\max_{n}\left|\mathcal{V}_n\right|\leq L\kappa$. Then,
\begin{equation}
\begin{array}{rl}
\rho_w\hspace{-2mm} &\displaystyle\leq \left|\left|\left(\widehat{\Sigma}_{Y_{w} Y_{w}}\right)^{\hspace{-1mm}-1} \widetilde{\Sigma}_{Y_{w} Y_{w}}\right|\right|_2 \\
&\displaystyle\leq \left|Y_w \right| \max_{i,i'}\left|\left[\left(\widehat{\Sigma}_{Y_{w} Y_{w}}\right)^{\hspace{-1mm}-1} \widetilde{\Sigma}_{Y_{w} Y_{w}}\right]_{ii'}\right|\\
&\displaystyle\leq KL^2\kappa\xi\varepsilon\ .
\end{array}
\label{rhobd}
\end{equation}
The first two inequalities follow from standard properties of matrix norm \citep{Golub96,Stewart90}. The last inequality is due to (\ref{max4}).

The rest of this proof utilizes the following result of \cite{Ipsen03} that is revised to reflect our notations:
\begin{thm}%[\citeauthor{Ipsen03}, \citeyear{Ipsen03}]
If 
%$\lambda_i>-1$ for $i=1,\ldots,|Y_{w}|$ and 
$|Y_{w}|\rho^2_w < 1$, then
%\begin{equation*}
$\log\hspace{-0.5mm}\left|\overline{\Sigma}_{Y_{w} Y_{w}}\right|\leq
\log\hspace{-0.5mm}\left|\widehat{\Sigma}_{Y_{w} Y_{w}}\right| \leq \log\hspace{-0.5mm}\left|\overline{\Sigma}_{Y_{w} Y_{w}}\right|-\log\hspace{-1mm}\left(1-|Y_{w}|\rho^2_w\right)
$
%\end{equation*}
for any joint walk $w$.
\label{bd1}
\end{thm} 
Using Theorem~\ref{bd1} followed by (\ref{rhobd}), 
\begin{equation}
\hspace{-6mm}
\begin{array}{rl}
\displaystyle\log\hspace{-0.5mm}\left|\widehat{\Sigma}_{Y_{w} Y_{w}}\right| - \log\hspace{-0.5mm}\left|\overline{\Sigma}_{Y_{w} Y_{w}}\right|
\hspace{-3mm} &\leq \displaystyle\log\hspace{-0.5mm}\frac{1}{1-|Y_{w}|\rho^2_w}\vspace{1mm}\\
\hspace{-3mm} &\leq \displaystyle\log\hspace{-0.5mm}\frac{1}{1\hspace{-0.5mm}-\hspace{-0.5mm}\left(K^{1.5}L^{2.5}\kappa\xi\varepsilon\right)^2}\vspace{-11mm}
\end{array}
\label{impt}
\end{equation}
for any joint walk $w$.
$$
\begin{array}{l}
\overline{\mathbb{H}}\hspace{-0.5mm}\left[Z_{Y_{w^{\ast}}}\right] - \overline{\mathbb{H}}\hspace{-0.5mm}\left[Z_{Y_{\widehat{w}}}\right]\vspace{1mm}\\
=\displaystyle\frac{1}{2}\left(
\left(\left|Y_{w^\ast}\right| - \left|Y_{\widehat{w}}\right|\right)\log\hspace{-0.5mm}\left(2\pi e\right) +
\log\hspace{-0.5mm}\left|\overline{\Sigma}_{Y_{w^\ast} Y_{w^\ast}}\right| - \log\hspace{-0.5mm}\left|\overline{\Sigma}_{Y_{\widehat{w}} Y_{\widehat{w}}}\right|\right)\vspace{1mm}\\
\leq \displaystyle\frac{1}{2}\left(
\left(\left|Y_{w^\ast}\right| - \left|Y_{\widehat{w}}\right|\right)\log\hspace{-0.5mm}\left(2\pi e\right) +
\log\hspace{-0.5mm}\left|\widehat{\Sigma}_{Y_{w^\ast} Y_{w^\ast}}\right| - \log\hspace{-0.5mm}\left|\overline{\Sigma}_{Y_{\widehat{w}} Y_{\widehat{w}}}\right|\right)\vspace{1mm}\\
\leq \displaystyle\frac{1}{2}\left(
\left(\left|Y_{\widehat{w}}\right| - \left|Y_{\widehat{w}}\right|\right)\log\hspace{-0.5mm}\left(2\pi e\right) +
\log\hspace{-0.5mm}\left|\widehat{\Sigma}_{Y_{\widehat{w}} Y_{\widehat{w}}}\right| - \log\hspace{-0.5mm}\left|\overline{\Sigma}_{Y_{\widehat{w}} Y_{\widehat{w}}}\right|\right)\vspace{1mm}\\
\leq \displaystyle\frac{1}{2}\log\frac{1}{1-\left(K^{1.5}L^{2.5}\kappa\xi\varepsilon\right)^2}\ .
\end{array}
$$
The first equality is due to (\ref{maxent21}).
% and $\left|Y_{w^\ast}\right| = \left|Y_{\widehat{w}}\right|$.
The first, second, and last inequalities follow from Theorem~\ref{bd1}, (\ref{bdiaga}),
and (\ref{impt}), respectively.
% and $\left|Y_{w^\ast}\right| = \left|Y_{\widehat{w}}\right|$
}{}
\end{document} 




uai2012.bbl
\begin{thebibliography}{}

\bibitem[Borg and Groenen(2005)Borg and Groenen]{borg05}
Borg, I. and Groenen, P. J.~F. (2005).
\newblock {\em Modern Multidimensional Scaling: Theory and Applications\/}.
\newblock Springer, NY.

\bibitem[Chen {\em et~al.\/}(2011)Chen, Rakha, and Sadek]{chen11}
Chen, H., Rakha, H.~A., and Sadek, S. (2011).
\newblock Real-time freeway traffic state prediction: A particle filter
  approach.
\newblock In {\em Proc. {IEEE ITSC}\/}, pages 626--631.

\bibitem[Chu {\em et~al.\/}(2007)Chu, Kim, Lin, Yu, Bradski, Ng, and
  Olukotun]{Chu07b}
Chu, C.-T., Kim, S.~K., Lin, Y.-A., Yu, Y., Bradski, G.~R., Ng, A.~Y., and
  Olukotun, K. (2007).
\newblock {Map-Reduce} for machine learning on multicore.
\newblock In B.~{Sch\"{o}lkopf}, J.~C. Platt, and T.~Hoffman, editors, {\em
  Advances in Neural Information Processing Systems 19\/}, pages 281--288,
  Cambridge, MA. MIT Press.

\bibitem[Chung {\em et~al.\/}(2004)Chung, Gupta, Burdick, and Murray]{Chung04}
Chung, T.~H., Gupta, V., Burdick, J.~W., and Murray, R.~M. (2004).
\newblock On a decentralized active sensing strategy using mobile sensor
  platforms in a network.
\newblock In {\em Proc. {CDC}\/}, pages 1914--1919.

\bibitem[Coates(2004)Coates]{Coates04}
Coates, M. (2004).
\newblock Distributed particle filters for sensor networks.
\newblock In {\em Proc. {IPSN}\/}, pages 99--107.

\bibitem[Cortes(2009)Cortes]{cortes09}
Cortes, J. (2009).
\newblock Distributed kriged {Kalman} filter for spatial estimation.
\newblock {\em {IEEE} Trans. Automat. Contr.\/}, {\bf 54}(12), 2816--2827.

\bibitem[Dolan {\em et~al.\/}(2009)Dolan, Podnar, Stancliff, Low, Elfes,
  Higinbotham, Hosler, Moisan, and Moisan]{LowSPIE09}
Dolan, J.~M., Podnar, G., Stancliff, S., Low, K.~H., Elfes, A., Higinbotham,
  J., Hosler, J.~C., Moisan, T.~A., and Moisan, J. (2009).
\newblock Cooperative aquatic sensing using the telesupervised adaptive ocean
  sensor fleet.
\newblock In {\em Proc. {SPIE} Conference on Remote Sensing of the Ocean, Sea
  Ice, and Large Water Regions\/}, volume 7473.

\bibitem[Golub and {Van Loan}(1996)Golub and {Van Loan}]{Golub96}
Golub, G.~H. and {Van Loan}, C.-F. (1996).
\newblock {\em Matrix Computations\/}.
\newblock Johns Hopkins Univ. Press, third edition.

\bibitem[Guestrin {\em et~al.\/}(2004)Guestrin, Bodik, Thibaus, Paskin, and
  Madden]{guestrin04}
Guestrin, C., Bodik, P., Thibaus, R., Paskin, M., and Madden, S. (2004).
\newblock Distributed regression: an efficient framework for modeling sensor
  network data.
\newblock In {\em Proc. {IPSN}\/}, pages 1--10.

\bibitem[Ipsen and Lee(2003)Ipsen and Lee]{Ipsen03}
Ipsen, I. C.~F. and Lee, D.~J. (2003).
\newblock Determinant approximations.
\newblock {Technical Report} CRSC-TR03-30, Center for Research in Scientific
  Computation, North Carolina State University, Raleigh, NC.

\bibitem[Kamarianakis and Prastacos(2003)Kamarianakis and
  Prastacos]{kamarianakis03}
Kamarianakis, Y. and Prastacos, P. (2003).
\newblock Forecasting traffic flow conditions in an urban network: Comparison
  of multivariate and univariate approaches.
\newblock {\em Transport. Res. Rec.\/}, {\bf 1857}, 74--84.

\bibitem[Krause {\em et~al.\/}(2008a)Krause, Horvitz, Kansal, and
  Zhao]{krause08}
Krause, A., Horvitz, E., Kansal, A., and Zhao, F. (2008a).
\newblock Toward community sensing.
\newblock In {\em Proc. {IPSN}\/}, pages 481--492.

\bibitem[Krause {\em et~al.\/}(2008b)Krause, Singh, and Guestrin]{Guestrin08}
Krause, A., Singh, A., and Guestrin, C. (2008b).
\newblock Near-optimal sensor placements in {Gaussian} processes: Theory,
  efficient algorithms and empirical studies.
\newblock {\em JMLR\/}, {\bf 9}, 235--284.

\bibitem[Lawrence {\em et~al.\/}(2003)Lawrence, Seeger, and
  Herbrich]{lawrence03}
Lawrence, N., Seeger, M., and Herbrich, R. (2003).
\newblock Fast sparse {Gaussian} process methods: The informative vector
  machine.
\newblock In S.~Becker, S.~Thrun, and K.~Obermayer, editors, {\em Advances in
  Neural Information Processing Systems 15\/}, pages 609--616, Cambridge, MA.
  MIT Press.

\bibitem[Low {\em et~al.\/}(2007)Low, Gordon, Dolan, and Khosla]{LowICRA07}
Low, K.~H., Gordon, G.~J., Dolan, J.~M., and Khosla, P. (2007).
\newblock Adaptive sampling for multi-robot wide-area exploration.
\newblock In {\em Proc. {IEEE ICRA}\/}, pages 755--760.

\bibitem[Low {\em et~al.\/}(2008)Low, Dolan, and Khosla]{LowAAMAS08}
Low, K.~H., Dolan, J.~M., and Khosla, P. (2008).
\newblock Adaptive multi-robot wide-area exploration and mapping.
\newblock In {\em Proc. {AAMAS}\/}, pages 23--30.

\bibitem[Low {\em et~al.\/}(2009a)Low, Dolan, and Khosla]{LowICAPS09}
Low, K.~H., Dolan, J.~M., and Khosla, P. (2009a).
\newblock Information-theoretic approach to efficient adaptive path planning
  for mobile robotic environmental sensing.
\newblock In {\em Proc. {ICAPS}\/}, pages 233--240.

\bibitem[Low {\em et~al.\/}(2009b)Low, Podnar, Stancliff, Dolan, and
  Elfes]{LowIPSN09}
Low, K.~H., Podnar, G., Stancliff, S., Dolan, J.~M., and Elfes, A. (2009b).
\newblock Robot boats as a mobile aquatic sensor network.
\newblock In {\em Proc. {IPSN-09} Workshop on Sensor Networks for Earth and
  Space Science Applications\/}.

\bibitem[Low {\em et~al.\/}(2011)Low, Dolan, and Khosla]{LowAAMAS11}
Low, K.~H., Dolan, J.~M., and Khosla, P. (2011).
\newblock Active {Markov} information-theoretic path planning for robotic
  environmental sensing.
\newblock In {\em Proc. {AAMAS}\/}, pages 753--760.

\bibitem[Low {\em et~al.\/}(2012)Low, Chen, Dolan, Chien, and
  Thompson]{LowAAMAS12}
Low, K.~H., Chen, J., Dolan, J.~M., Chien, S., and Thompson, D.~R. (2012).
\newblock Decentralized active robotic exploration and mapping for
  probabilistic field classification in environmental sensing.
\newblock In {\em Proc. {AAMAS}\/}, pages 105--112.

\bibitem[Min and Wynter(2011)Min and Wynter]{min11}
Min, W. and Wynter, L. (2011).
\newblock Real-time road traffic prediction with spatio-temporal correlations.
\newblock {\em Transport. Res. C-Emer.\/}, {\bf 19}(4), 606--616.

\bibitem[Neumann {\em et~al.\/}(2009)Neumann, Kersting, Xu, and
  Schulz]{Kersting09}
Neumann, M., Kersting, K., Xu, Z., and Schulz, D. (2009).
\newblock Stacked {Gaussian} process learning.
\newblock In {\em Proc. {ICDM}\/}, pages 387--396.

\bibitem[Oh {\em et~al.\/}(2010)Oh, Xu, and Choi]{oh10}
Oh, S., Xu, Y., and Choi, J. (2010).
\newblock Explorative navigation of mobile sensor networks using sparse
  {Gaussian} processes.
\newblock In {\em Proc. CDC\/}, pages 3851--3856.

\bibitem[{Olfati-Saber}(2005){Olfati-Saber}]{saber05a}
{Olfati-Saber}, R. (2005).
\newblock Distributed {Kalman} filter with embedded consensus filters.
\newblock In {\em Proc. {CDC}\/}, pages 8179--8184.

\bibitem[Paskin and Guestrin(2004)Paskin and Guestrin]{Guestrin04b}
Paskin, M.~A. and Guestrin, C. (2004).
\newblock Robust probabilistic inference in distributed systems.
\newblock In {\em Proc. {UAI}\/}, pages 436--445.

\bibitem[Podnar {\em et~al.\/}(2010)Podnar, Dolan, Low, and
  Elfes]{LowAeroconf10}
Podnar, G., Dolan, J.~M., Low, K.~H., and Elfes, A. (2010).
\newblock Telesupervised remote surface water quality sensing.
\newblock In {\em Proc. {IEEE} Aerospace Conference\/}.

\bibitem[{Qui{\~{n}}onero-Candela} and Rasmussen(2005){Qui{\~{n}}onero-Candela}
  and Rasmussen]{candela05}
{Qui{\~{n}}onero-Candela}, J. and Rasmussen, C.~E. (2005).
\newblock A unifying view of sparse approximate {Gaussian} process regression.
\newblock {\em JMLR\/}, {\bf 6}, 1939--1959.

\bibitem[Rasmussen and Williams(2006)Rasmussen and Williams]{rasmussen06}
Rasmussen, C.~E. and Williams, C. K.~I. (2006).
\newblock {\em Gaussian Processes for Machine Learning\/}.
\newblock MIT Press, Cambridge, MA.

\bibitem[Rosencrantz {\em et~al.\/}(2003)Rosencrantz, Gordon, and
  Thrun]{Thrun03}
Rosencrantz, M., Gordon, G., and Thrun, S. (2003).
\newblock Decentralized sensor fusion with distributed particle filters.
\newblock In {\em Proc. {UAI}\/}, pages 493--500.

\bibitem[Sch{\"{o}}lkopf and Smola(2002)Sch{\"{o}}lkopf and Smola]{scholkopf02}
Sch{\"{o}}lkopf, B. and Smola, A.~J. (2002).
\newblock {\em Learning with Kernels: Support Vector Machines, Regularization,
  Optimization, and Beyond\/}.
\newblock MIT Press, Cambridge, MA, first edition.

\bibitem[Schrank {\em et~al.\/}(2011)Schrank, Lomax, and Eisele]{Schrank11}
Schrank, D., Lomax, T., and Eisele, B. (2011).
\newblock {\em {TTI}'s 2011 Urban Mobility Report\/}.
\newblock {Texas Transportation Institute, Texas A\&M University}.

\bibitem[Seeger and Williams(2003)Seeger and Williams]{seeger03}
Seeger, M. and Williams, C. (2003).
\newblock Fast forward selection to speed up sparse {G}aussian process
  regression.
\newblock In C.~M. Bishop and B.~J. Frey, editors, {\em Proc. {AISTATS}\/}.

\bibitem[Srinivasan and Jovanis(1996)Srinivasan and Jovanis]{Srinivasan96}
Srinivasan, K.~K. and Jovanis, P.~P. (1996).
\newblock Determination of number of probe vehicle required for reliable travel
  time measurement in urban network.
\newblock {\em Transport. Res. Rec.\/}, {\bf 1537}, 15--22.

\bibitem[Stewart and Sun(1990)Stewart and Sun]{Stewart90}
Stewart, G.~W. and Sun, J.-G. (1990).
\newblock {\em Matrix Perturbation Theory\/}.
\newblock Academic Press.

\bibitem[Stranders {\em et~al.\/}(2009)Stranders, Farinelli, Rogers, and
  Jennings]{stranders09}
Stranders, R., Farinelli, A., Rogers, A., and Jennings, N.~R. (2009).
\newblock Decentralised coordination of mobile sensors using the max-sum
  algorithm.
\newblock In {\em Proc. IJCAI\/}, pages 299--304.

\bibitem[Sukkarieh {\em et~al.\/}(2003)Sukkarieh, Nettleton, Kim, Ridley,
  Goktogan, and {Durrant-Whyte}]{Durrant-Whyte03b}
Sukkarieh, S., Nettleton, E., Kim, J., Ridley, M., Goktogan, A., and
  {Durrant-Whyte}, H. (2003).
\newblock The {ANSER} project: Data fusion across multiple uninhabited air
  vehicles.
\newblock {\em IJRR\/}, {\bf 22}(7-8), 505--539.

\bibitem[Turner {\em et~al.\/}(1998)Turner, Eisele, Benz, and
  Holdener]{turner98}
Turner, S.~M., Eisele, W.~L., Benz, R.~J., and Holdener, D.~J. (1998).
\newblock Travel time data collection handbook.
\newblock {Technical Report} FHWA-PL-98-035, Federal Highway Administration,
  Office of Highway Information Management, Washington, DC.

\bibitem[Wang and Papageorgiou(2005)Wang and Papageorgiou]{wang05}
Wang, Y. and Papageorgiou, M. (2005).
\newblock Real-time freeway traffic state estimation based on extended {Kalman}
  filter: a general approach.
\newblock {\em Transport. Res. B-Meth.\/}, {\bf 39}(2), 141--167.

\bibitem[Work {\em et~al.\/}(2010)Work, Blandin, Tossavainen, Piccoli, and
  Bayen]{Bayen10a}
Work, D.~B., Blandin, S., Tossavainen, O., Piccoli, B., and Bayen, A. (2010).
\newblock A traffic model for velocity data assimilation.
\newblock {\em AMRX\/}, {\bf 2010}(1), 1--35.

\end{thebibliography}
