A Differentially Private Algorithm for Range Queries on Trajectories

TL;DR

This paper introduces a new differentially private algorithm for range queries on trajectory data that adaptively adds noise, achieving lower error and higher utility than existing methods by considering data and query awareness.

Contribution

It presents the first data- and query-aware differential privacy method for trajectory range queries that adaptively partitions data to improve accuracy.

Findings

Achieves significantly lower error compared to state-of-the-art methods.

Effectively maintains trajectory data consistency while ensuring privacy.

Demonstrates high accuracy and efficiency on real and synthetic datasets.

Abstract

We propose a novel algorithm to ensure $\epsilon$-differential privacy for answering range queries on trajectory data. In order to guarantee privacy, differential privacy mechanisms add noise to either data or query, thus introducing errors to queries made and potentially decreasing the utility of information. In contrast to the state-of-the-art, our method achieves significantly lower error as it is the first data- and query-aware approach for such queries. The key challenge for answering range queries on trajectory data privately is to ensure an accurate count. Simply representing a trajectory as a set instead of \emph{sequence} of points will generally lead to highly inaccurate query answers as it ignores the sequential dependency of location points in trajectories, i.e., will violate the consistency of trajectory data. Furthermore, trajectories are generally unevenly distributed across a city and adding noise uniformly will generally lead to a poor utility. To achieve differential privacy, our algorithm adaptively adds noise to the input data according to the given query set. It first privately partitions the data space into uniform regions and computes the traffic density of each region. The regions and their densities, in addition to the given query set, are then used to estimate the distribution of trajectories over the queried space, which ensures high accuracy for the given query set. We show the accuracy and efficiency of our algorithm using extensive empirical evaluations on real and synthetic data sets.

Tables

Table 1. TABLE I : Summary of notation.

$D$	Trajectory data set, $D\subseteq 𝒯$ and $\left|D\right|=n$
$H$	Spatial histogram related to $D$, $H=\{F,E\}$, $H\subseteq \mathscr{H}$
$E$	The horizontal and vertical edges in $H$, $E=\{E_v,E_h\}$
$j$	Index of a cell in histogram, $j=(row,col)$
$H\left[j\right]$	A cell in $H$ with a face $f\in F$ and four edges $e\in E$
$f_j$	The face component of cell $H\left[j\right]$
$e_{ij}$	The incident edge of two adjacent cells $H\left[i\right]$ and $H\left[j\right]$
$H^{\prime }$	Synthetic spatial histogram
$H^{\prime \prime }$	Intermediate synthetic histogram in the synthesizing process
$T$	Total number of iteration in the synthesizing process
$Q$	Set of queries $q\in Q$ that $q=\{q_F,q_E\}$
$q\left[j\right]$	The $j$th element of query, $q\left[j\right]\in \{0,1\}$
$q\left(H\right)$	Computing the query $q$ on histogram $H$
$P$	Set of partitions $p_i\in P$
$p_i\left[j\right]$	The $j$th cell in the region $p_i$
$B$	Set of partition densities $b_i\in B$, $b_i={\sum }_{j}c\left(p_i\left[j\right]\right)/\left|D\right|$
$\delta$	The threshold of uniformity cost in partitioning (Algorithm 1)
$wer{r}_i$	Error of selected query at iteration $i$ (Algorithm 2)
$c\left(x\right)$	The count of $x$
$\left|x\right|$	Total number of elements/cells in $x$

Table 2. TABLE II : DQAM versus existing works in terms of technical and output features.

Mechanism	Data-aware	Query-aware	Greedily Consistent	Optimally Consistent
LM		✓
MWEM		✓
DAWA	✓	✓
PriSH		✓	✓
DQAM	✓	✓	✓	✓

Table 3. TABLE III : Standard deviation of DQAM error and its competing algorithms for each ϵ italic-ϵ \epsilon on the sample of Porto-Taxi with 1000 1000 1000 trajectories.

$ϵ$	LM	MWEM	DAWA	PriSH	DQAM
$0.01$	500.86	5.64	3.65	2.50	1.73
$0.05$	500.15	4.48	2.00	2.81	1.54
$0.1$	429.71	3.37	1.24	1.67	0.53
$0.5$	285.44	3.10	0.86	0.39	0.27

Equations

$$q(H)=q_F(F)-q_E(E),$$

$$c(H[j]{)}_{j\in p}=\frac{{\sum }_{i\in p}c(H[i])}{|p|},$$

$$cost(p)=\sum _{j\in p}\mid c(H[j])-\frac{{\sum }_{i\in p}c(H[i])}{|p|}\mid .$$

$$b_i=\frac{{\sum }_{j\in p_i}c(H[j])}{n}.$$

$$cost(p)=cost(p)+\nu .$$

$$cos{t}^{\prime }(p)$$

$$+$$

$$\le$$

$$=$$

$$\Delta cost=cos{t}^{\prime }(p)-cost(p)\le 2k.$$

$$density(p)=\frac{{\sum }_{j\in p}c(H[j])}{n}+{\nu }^{\prime }.$$

$$s(q_i,H)=|q_i(H)-q_i(H_i^{\prime })|.$$

$$H_{i+1}^{\prime }[j]\propto H_i^{\prime }[j].exp(q_i[j].wer{r}_i.b_p/2n).$$

$$wer{r}_i=(q(H)+{\nu }_i)-q(H_i^{\prime }).$$

$$s_i\propto exp(q_i[j].wer{r}_i).$$

$$s_i\propto exp(q_i[j].wer{r}_i.b_p),$$

$$H_{i+1}^{\prime }[j]\propto H_i^{\prime }[j].exp(q_i[j].wer{r}_i.b_p/2n),$$

$$c(e_{ij})\le \min \{c(f_i),c(f_j)\},\forall e_{ij}\in E^{\prime \prime },\text{and}\forall f_i,f_j\in F^{\prime \prime }.$$

$$\underset{H^{\prime \prime }}{\mathrm{minimize}}$$

$$c(H^{\prime \prime }[i])\ge 0,\forall i\in \{1,...,|H^{\prime \prime }|\}$$

$$c(e_{ij})\le \min \{c(f_i),c(f_j)\},\forall e_{ij}\in E^{\prime \prime },\forall f_i,f_j\in F^{\prime \prime }.$$

$$\underset{q\in Q}{\max }|q(H)-q(H^{\prime })|\le 2n\sqrt{\frac{log|H|}{T}}+\frac{10Tlog|Q|}{ϵ}$$

$$RE(H||H^{\prime })=\sum _{j\in H}c(H[j])log(\frac{c(H[j])}{c(H^{\prime }[j])})/n,$$

$$P(\mathcal{A}\in L_{\rho })$$

$$=$$

$$\ge$$

$$=$$

$$\alpha =\underset{i\in T}{\max }\underset{x\in H^{\prime }}{\max }{\alpha }_x.$$

$${\psi }_{i-1}-{\psi }_i=\sum _{x\in H}c(H[x])log(\frac{c(H_i^{\prime }[x])}{c(H_{i-1}^{\prime }[x])})/n,$$

$${\psi }_{i-1}-{\psi }_i$$

$${\beta }_i$$

$${\psi }_{i-1}-{\psi }_i$$

$${\psi }_{i-1}-{\psi }_i\ge (\frac{q_i(H_{i-1}^{\prime })-q_i(H)}{2n}{)}^2-(\frac{m_i-q_i(H)}{2n}{)}^2,$$

Full text

I Introduction

The popularity of sensor-enabled devices (e.g., wearables and smartphones) has significantly advanced the capability of businesses to collect and analyze people’s trajectories. Studying large-scale data sets and analyzing the movement patterns of individuals provide crucial insight for many applications (e.g., traffic management, route planning, urban planning, crime detection). Such applications critically rely on estimating the number of trajectories in an area. For example, in urban planning, computing the number of pedestrians and the flow of their movements provide significant information for the placement of public spaces (e.g., local parks, street spaces, plazas), pedestrian and bicycle paths, and public transport stations. A fundamental query type in studying the movement patterns is range query on trajectories [12, 14] which counts the number of distinct trajectories intersecting the query area (2D space). However, computing the cost of this query type on raw trajectories is linear in the number of trajectories and the average length of the trajectories. To avoid such high computational cost spatial histogram (see Section II) is commonly used to compute range queries efficiently. The cost of range queries using spatial histogram is linear in the size of the query area. A spatial histogram considers the data space as a grid with fixed cell size and aggregates the count of trajectories overlapping a cell or an incident edge of two cells in the grid. In other words, a spatial histogram maps the trajectory data into a set of counts representing the distribution of trajectories throughout the space.

Even though the value of spatial histograms in studying movement patterns is realized by companies, releasing such data to third parties is still a significant issue due to privacy concerns. Studies show that trajectories exhibit a high risk of being identified that aggregation or coarsening techniques as done in a spatial histogram, barely reduce their uniqueness [25, 4]. The key reason is that trajectories are time-stamped sequences of locations with a start and end points. For example, the start of a trajectory is often a home location. Thus, publishing a spatial histogram raises privacy concerns due to the uniqueness of people’s trajectories, particularly in sparse regions with few trajectory counts.

A principled approach is needed to generate a synthetic spatial histogram with strong privacy guarantees of differential privacy [5] while ensuring high utility in practice. The utility refers to the accuracy in answering range queries. The approach can make use of a set of range queries to learn the distribution of data (called query-awareness). However, it should ensure the synthesized spatial histogram maintains high accuracy in answering any range queries on trajectories. A naive approach adds an equal noise to the cell counts in the histogram to ensure privacy. Since a trajectory may intersect multiple cells, the cells are sequentially dependent. Thus, naively scaling noise to satisfy differential privacy can violate the sequentiality between the cells, which is crucial in representing trajectories. Furthermore, the sensitivity of cell counts to the added noise varies depending on the density of trajectories in different regions of data space. A large region of space might be sparse with small cell counts which makes the cell counts in the region highly sensitive to the added noise. In contrast, a region might be small but very dense such that cell counts are large and highly noise resistance. The proposed approach should be able to ensure maintaining the data properties (called data-awareness) in the synthesized data.

In this paper, we present Data- and Query-Aware Mechanism (DQAM), a mechanism that synthesizes spatial histograms. As [6], DQAM takes a spatial histogram and a set of range queries as input and uses the correlation between the queries (i.e., is query-aware) to synthesize the histogram. In contrast to [6], DQAM identifies the density of different regions in data and make use of the density in scaling the added noise to the generated synthetic histogram (i.e., is query-aware). Also, DQAM generates the optimal synthetic histogram in each iteration of the synthesizing process that satisfies the sequentiality of cells. The proposed mechanism in [6] applies a greedy strategy to satisfy the sequentiality requirement among histogram cells, which may lead to overcorrection. While the work in [6] successfully deals with trajectories uniformly distributed in the data space, it can result in considerable loss of utility for trajectory data sets with a non-uniform distribution in different regions. DQAM is designed to deal with such realistic trajectory distributions by capturing the density of trajectories in different regions and using the density as a weighting measure: each region is assigned a weight relative to its density. To guarantee differential privacy, DQAM uses the well established Laplace mechanism [5] and Exponential mechanism [15] for adding noise in all parts of its computations.

Contributions. First, we propose a differentially private Data- and Query-Aware Mechanism (DQAM), a mechanism that publishes a spatial histogram to answer range queries efficiently. To the best of our knowledge, DQAM is the first data- and query-aware mechanism for range queries on sequential spatial data that reduces the error by a factor up to 7.4 compared to the current approaches on trajectory data sets. Second, we design an efficient data-aware algorithm with $nlog(n)$ time complexity in the number of cells in the spatial histogram that partitions the histogram into uniform regions. Third, we present a query-aware and differentially private strategy that captures trajectory sequentiality and synthesizes accurate output using the given queries and partitions.

Mechanism Overview

Differential privacy assumes a privacy budget for an algorithm that is described by a parameter $\epsilon$. DQAM is a two-stage $\epsilon$-differentially private mechanism for answering range queries on trajectories (see Fig. 1). The mechanism takes as input a set of range queries $Q$, and a spatial histogram $H$. The histogram $H$ is a representation of trajectory counts throughout the data space. The query set $Q$ is randomly generated and is large to cover the entire space of $H$. The output of DQAM is $H^{\prime}$, a private estimate of $H$, which preserves the distribution of original histogram. In first stage, the histogram is partitioned into disjoint and uniform regions and the trajectory density of each region is computed. Intuitively, partitioning the histogram space into uniform regions improves the signal to noise ratio, which in turn improves the accuracy of the estimation process in stage 2. Instead of estimating the density of every single cell, the algorithm estimates the density of a region and assumes the same density for all cells in the region. Additionally, maintaining the density of a uniform region is significantly easier than a non-uniform region. Since the trajectory data sets are usually sparse, there are potentially many regions with (semi-)uniform distribution in the data. Utilizing this property, the estimation algorithm considers a uniform region as a single entity and uses its density in the private estimation process. Furthermore, the requirement that generated partitions should be disjoint, ensures that each region can be treated independently.

Any algorithm that efficiently generates disjoint and uniform partitions can be used in this stage of DQAM. One approach is to generate all possible partitions and compute the uniformity cost of each partition. The partitioning algorithm should find a set of these partitions, which are disjoint, cover the entire space, and have the minimum uniformity cost. However, there is no closed-form expression for computing all possible partitions in a $m\times n$ space and the computational complexity grows with the order of $O(2^{m\times n})$. Furthermore, finding the set with the minimum cost is computationally expensive. Another approach is using space partitioning techniques such as quadtree, kd-tree and R-tree. However, the R-tree and kd-tree cannot adapt well to identify disjoint and uniform regions of trajectories. An R-tree cannot cover areas with no trajectory, and the output partitions may be overlapped. A kd-tree requires an ordering among trajectories to do partitioning while neither trajectories nor their counts in the histogram have any order. In addition, the partitioning in kd-tree is based on the density of regions while the outputs might not have a uniform distribution. In contrast, a quadtree decomposes the space into disjoint regions. By defining a partitioning cost to measure the uniformity of regions, a quadtree can partition the space into uniform regions. Whilst, the generated partitions may not be optimal concerning the uniformity cost, the quadtree is computationally efficient, and the final partitioning provides an accurate estimation of uniform regions. However, the regions in a quadtree branch may overlap which could be used to refine parts of a user’s trajectory. We need to consider this correlation in the variance of added noise to ensure differential privacy.

Note that regardless of the approach used in stage one, only the outputs (i.e., final partitioning) are used for the next stage of DQAM. Thus, we need to ensure that the final partitions are differentially private as the rest of the generated partitions during the process will be removed. For example, in a constructed quadtree, we need to ensure the privacy of leaves.

In the second stage, DQAM uses the computed final partitions and densities, in addition to the given query set $Q$, to estimate the distribution of original histogram $H$. The output is an $\epsilon$-differentially private histogram $H^{\prime}$ that is generated using the data and query properties, i.e., DQAM is data- and query-aware. Since the two stages iteratively access the original histogram, we split the $\epsilon$ budget. We will show that we need in total four privacy budgets $\epsilon_{i},i\in\{1,\ldots,4\}$ with $\sum_{i=1}^{4}\epsilon_{i}=\epsilon$. We conducted a pre-study that investigated the impact of different weightings but found no significant impact on the accuracy/utility and hence, set $\epsilon_{i}=\frac{\epsilon}{4}$.

Step 1: Private Partitioning. We develop an efficient algorithm based on the two approaches described above that defines the uniformity as a cost function and partitions the histogram into uniform regions. Intuitively, in a uniform region, all values can be manipulated/treated equally, which significantly improves the accuracy of the synthesizing process (see Section III-A). It generates a list of uniform regions $P$ (i.e., the quadtree leaves) representing disjoint partitions of the histogram. The outputs of this stage are the partitions ($P$) and their densities ($B$), i.e., the number of trajectories in a region relative to the total number of trajectories in the entire histogram. Since the partitioning process and computing densities need interacting with the original histogram, we use the Laplace mechanism [5] to add noise to them.

Step 2: Private Synthesizing Process. In this step, our algorithm takes the partitions $P$, densities $B$ and the given range queries $Q$ to privately learn the distribution of the original spatial histogram $H$ and generate an $\epsilon$-differentially private histogram $H^{\prime}$. The learning process is iterative and starts with a uniform trajectory distribution. In each iteration, one query $q\in Q$ is selected. Based on the query error on $H^{\prime}$, the estimation is updated by rescaling up/down the values in $H^{\prime}$. The values in a region $p\in P$ are rescaled if $p$ overlaps with $q$. The density of region $p$ magnifies/reduces the update scale for that region. Previous work on estimating histogram values while ensuring local sequentiality [6] used a heuristic strategy to preserve the dependency. Our paper presents an algorithm to compute the optimal estimation for $H^{\prime}$ while ensuring the dependency. For the query selection and adding noise to the query answer, the exponential and Laplace mechanisms are applied, respectively.

The following example is a sample execution of DQAM.

Example 1. An overview of DQAM is presented in Fig. 1. The inputs are a histogram $H$ with $4$ trajectories and the set of range queries $Q$ shown graphically as rectangular areas over the histogram. A possible output of partitioning can be $P=\{p_{1},...,p_{7}\}$. This need not to be optimal as defined in Section III because the partitioning is randomized. The density of partitions are $B=\{b_{1}=2/4=0.5,b_{2}=1/4=0.25,b_{3}=1/4=0.25,b_{4}=1/4=0.25,b_{5}=3/4=0.75,b_{6}=1/4=0.25,b_{7}=1/4=0.25\}$. The $P$ and $B$ in addition to $Q$ are then used in step 2 to generate $H^{\prime}$.

II Background

II-A Spatial Histogram and Range Queries

Consider a range query with a convex area over a trajectory data set $D$ which counts the distinct number of trajectories intersecting the area. Let $D$ be a subset of size $n$ of the trajectory universe $\mathcal{T}$ of trajectories without loops. Each trajectory $t\in D$ has a domain equal to the data space $\mathcal{S}=\{s_{1},s_{2},...,s_{N}\}$ where $s_{i}\in\mathcal{S}$ is a location given by longitude and latitude. To express the query, we first map $D$ into a grid-based histogram called spatial histogram (see Fig. 2).

A spatial histogram $H=\{F,E_{v},E_{h}\}$ with $r$ rows and $c$ columns, consist of three components: (1) faces $F\in\mathbb{N}_{0}^{r\times c}$, (2) vertical edges $E_{v}\in\mathbb{N}_{0}^{r\times(c-1)}$ and (3) horizontal edges $E_{h}\in\mathbb{N}_{0}^{(r-1)\times c}$. In particular, $H$ maps $D$ into three matrices describing the overlap of trajectories with the histogram cells, i.e., $H\mathrel{\mathop{\mathchar 58\relax}}D\rightarrow\mathcal{H}$, where $\mathcal{H}=\{\mathbb{N}_{0}^{r\times c},\mathbb{N}_{0}^{r\times(c-1)},\mathbb{N}_{0}^{(r-1)\times c}\}$ is the histogram domain. The edge components (small boxes in Fig. 2) distinguish a spatial histogram from simple histograms by recording the count of trajectories crossing the edge between two adjacent cells. The edge components preserve the sequential dependency of locations in a trajectory referred to as local sequentiality. We use $E=\{E_{v},E_{h}\}$ from now on to refer to both the edge components. To make the presentation simple, we use $e\in E$ to denote an edge $e$ in which $e\in E_{v}$ or $e\in E_{h}$. Without loss of generality we assume each histogram cell $H[j]$ in index $j=(row,col)$ is a square of fixed side length and all cells have the same size.

Spatial histograms are known as data structures to efficiently store trajectories and compute range queries [23, 24]. The shaded areas $q_{1}$ and $q_{2}$ in Fig. 2 show two range queries. We assume a range query specifies a rectangular region of $r^{\prime}\geq r$ cells height and $c^{\prime}\geq c$ cells width, where the axes are aligned with the cells’ boundaries. A range query can be answered by subtracting the sum of edge values from the sum of face values, which are both inside the query area. For instance, the total number of trajectories in $q_{2}$ is $2-1=1$. Note the border edges are not counted in computing the query answer because spatial histograms are constructed in such a way that trajectories can only start and end in faces.

Definition II.1**.**

A range query $q=\{q_{F},q_{E}\mid q_{F}=$ set of faces inside the query rectangle q and $q_{E}=$ set of edges inside the query rectangle $q\}$ is as a function that takes a histogram $H=\{F,E\}$ as input and maps it to a non-negative integer values, i.e., $q\mathrel{\mathop{\mathchar 58\relax}}\mathcal{H}\rightarrow\mathbb{N}_{0}$, by counting the distinct trajectories intersecting the query area:

[TABLE]

where $q_{F}(F)=q_{F}\diamond F$ and $q_{E}(E)=q_{E}\diamond E$, and “$\diamond$” is the dot product defined over matrices (i.e., the sum of element-wise multiplication of the matrices. This can also be written as trace of $(q_{F}*F^{\intercal})$).

Note that $q_{F}$ and $q_{E}$ are binary matrices indicating the query area. Our definition would translate to the general definition of a range query if we used the lower and upper column and row indices. A trajectory may have multiple intersections with the query area (e.g., $q_{1}$ in Fig. 2). Computing the range query for such cases requires modified versions of spatial histograms [22, 24]. For the sake of simplicity, we use the basic histogram defined above, but the results can be generalized. The notations of the paper are summarized in Table I.

II-B Differential Privacy

According to differential privacy [5], the difference between the probability of mechanism outputs, in the presence/absence of a record in the input data set, is bounded by a privacy budget $\epsilon$. For a trajectory data set $D$, the set $nbr(D)$ denotes data sets that are different from $D$ in at most one trajectory (by changing, adding or removing a trajectory), i.e., nbr(D)=\{D^{*}\big{|}|D-D^{*}|=1\}. We refer to the spatial histogram $H$ instead of $D$, since $H$ represents a mapping of the trajectory data set $D$ and changing a trajectory in $D$ directly affects the counts in $H$.

Definition II.2**.**

Given a histogram $H\subseteq\mathcal{H}$ and any histogram $H^{*}\in nbr(H)$, a randomized mechanism $\mathcal{M}$ is $\epsilon$-differentially private, $\epsilon>0$, and any output $O\subseteq Range(\mathcal{M})$ if $Pr[\mathcal{M}(H)\in O]\leq exp(\epsilon).Pr[\mathcal{M}(H^{*})\in O]$.

Differential privacy has two important composition properties [16]. Let $\mathcal{M}_{1},\mathcal{M}_{2},...\mathcal{M}_{n}$ be $n$ mechanisms each satisfying $\epsilon_{i}$-differential privacy. When $\mathcal{M}_{1},\mathcal{M}_{2},...\mathcal{M}_{m}$ are executed sequentially on a histogram, it is called sequential composition and the output is $(\sum\epsilon_{i})$-differentially private. When the histogram is partitioned into $n$ disjoint subsets and $\mathcal{M}_{i}$ is applied on subset $i$, it is called parallel composition and the output is $\{\max{\epsilon_{i}}\}$-differentially private. Since in a spatial histogram a trajectory is represented as a sequence of cells, each count can be treated separately, i.e., parallel composition. However, a naive approach may violate the local sequentiality of cells, which significantly reduces the utility of the output. Our proposed mechanism guarantees the local sequentiality in the released histogram and thus achieves high utility.

For a range query $q$ differential privacy can be achieved by adding a properly scaled random noise to the output. The noise scale depends on $q$’s sensitivity [5], which captures the maximum difference in query answers between any two neighboring histograms:

Definition II.3**.**

The sensitivity of $q\mathrel{\mathop{\mathchar 58\relax}}\mathcal{H}\rightarrow\mathbb{N}_{0}$, denoted by $\Delta q$, is $\Delta q=\max\limits_{H,H^{*}\in nbr(H)}||q(H)-q(H^{*})||_{1}$.

Since adding (removing) one trajectory to (from) a spatial histogram changes the answer of a range query at most by one, the sensitivity of a range query is $1$. The Laplace mechanism (LM) [5] achieves differential privacy by adding a Laplace noise to the query’s answer. We use Lap($\sigma$) to denote the Laplace probability distribution with mean [math] and scale $\sigma$.

Definition II.4**.**

Given a histogram $H\subseteq\mathcal{H}$ and a range query $q\mathrel{\mathop{\mathchar 58\relax}}\mathcal{H}\rightarrow\mathbb{N}_{0}$, the Laplace mechanism is defined as $LM(H,q)=q(H)+Lap(\Delta q/\epsilon)^{|H|}$.

We apply the Laplace mechanism in the first step of DQAM where the cost and density of a partition are functions of a range query. The second step, however, is a non-numerical function that selects a query $q$ from $Q=\{q_{i}|i\in\{1,2,...,m\}\}$. To achieve differential privacy, the Exponential Mechanism (EM) [15] randomly selects a query with the highest utility on the estimated histogram $H^{\prime}$. In the Exponential mechanism, a random noise, which is scaled to the sensitivity of the utility function, is added to the utility scores [15]. The sensitivity bounds the difference between the probability of selecting $q$ on any two neighbouring histograms.

Definition II.5**.**

The sensitivity of a utility function $u\mathrel{\mathop{\mathchar 58\relax}}\mathcal{H}\times Q\rightarrow\mathbb{R}$ to score each query $q\in Q$, denoted by $\Delta u\in[0,1]$, is $\Delta u=\max_{q\in Q}\max_{H^{\prime},H^{\prime*}\in nbr(H^{\prime})}||u(H^{\prime},q)-u(H^{\prime*},q)||_{1}$.

Definition II.6**.**

Given a histogram $H\subseteq\mathcal{H}$, a query set $Q$ and a utility function $u\mathrel{\mathop{\mathchar 58\relax}}\mathcal{H}\times Q\rightarrow\mathbb{R}$, the exponential mechanism $EM(H,Q,u)$ outputs a query $q\in Q$ proportional to $exp(\frac{\epsilon u(H,q)}{2\Delta u})$.

By bounding the sensitivity ($\Delta q$ in Laplace mechanism or $\Delta u$ inExponential mechanism), the privacy budget $\epsilon$ identifies the scale of injected noise.

III Step 1: Private Partitioning

The first step of DQAM partitions space into uniform regions $P$ and computes their densities $B$. This stage is independent of the query set and finds $P$ and $B$ as two parameters describing data properties. These parameters are later used in the synthesizing process (step two) to estimate $H^{\prime}$ as close as possible to $H$ without violating $\epsilon$-differential privacy. $\epsilon_{1}$ is used for computing the cost of a partition and $\epsilon_{2}$ for injecting Laplace noise to the partitions densities, where $\epsilon_{1}+\epsilon_{2}=\epsilon/2$. Algorithm 1 describes our private partitioning and is discussed below.

III-A Quadtree search

As discussed earlier, the goal is to find partitions in the histogram that are disjoint and uniform. These requirements are needed to minimize the error in the update step (see s2 in Algorithm 2) where we scale all counts in a region with the same magnitude. To achieve this goal, we propose an approach which is a combination of the two approaches described earlier in Section I. In particular, instead of generating all possible partitions in the histogram, we approximate them with all partitions that could be generated if the histogram was decomposed using a quadtree. Whilst it might come at the expense of a small cost of utility, it leads to a highly efficient algorithm: (1) Given the partitions in different levels of a quadtree, we can compute the uniformity cost of each partition and add noise to ensure differential privacy; (2) Next, a search procedure expands a quadtree structure can efficiently evaluate all possible combinations and find the least cost set of partitions that are disjoint and cover the entire space of histogram.

The procedure shown in Algorithm 1 combines the two stages. It takes $H$ as input, and returns the leaves of a (potentially unbalanced) quadtree $P=\{p_{1},...,p_{k}\}$ with their densities $B=\{b_{1},...,b_{k}\}$. Using the full histogram $H$ as the root node, the algorithm recursively splits $H$ into $4$ regions with a size of power of two using the split() function. The noisy cost of each region in the current level (pcost) is compared to the average of noisy cost of its four children (chcost). If the difference between chcost and pcost is greater than a threshold $\delta$, the partitioning proceeds to the next level. Otherwise, the node is labeled as “leaf”. Since computing the costs require interacting with the original histogram, a random value $\nu\sim Lap(\frac{\Delta cost}{\epsilon_{1}})$ is added to each region’s cost. The algorithm terminates when all regions are labeled as “leaf”, or the size of each region is $1$ (one cell). Intuitively, when the node has a good level of uniformity, the children would not significantly improve the uniformity cost.

We define the partitioning cost function to measure the uniformity of a region. Intuitively, a high trajectory density in a part of region incurs a high cost and triggers the algorithm to split the region into smaller ones. Let $j$ denote the cell index in the histogram $H$ and region $p$ be a subset of cell indexes of $H$. We say a region $p$ is uniform when

[TABLE]

where $|\cdot|$ and $c(\cdot)$ represent the set size and the cell count, respectively. Fig. 3 shows examples of uniform and non-uniform regions. With this observation, the cost of a region $p$, denoted by cost(p), is the amount that $p$ deviates from being perfectly uniform:

[TABLE]

III-B Computing densities

After privately identifying uniform regions $P=\{p_{1},...,p_{k}\}$, we need to compute the densities of regions $B=\{b_{1},...,b_{k}\}$. The density of region $p_{i}\in P$, named $b_{i}$ defines as the number of trajectories in $p_{i}$ relative to the total number trajectories in the histogram. Given $n$, the total number of trajectories in histogram, the density $b_{i}$ of region $p_{i}$ is:

[TABLE]

The function $density(H)$ in Algorithm 1 denotes computing the density of regions in the histogram $H$.

In the second step of DQAM, a region specifies the area to be updated, and its density describes the update magnitude.

III-C Formal guarantees

We analyze Algorithm 1 for three key aspects: privacy, accuracy and efficiency.

Privacy. We show that step 1 is $\frac{\epsilon}{2}$-differentially private.

Theorem III.1**.**

Algorithm 1 satisfies $\frac{\epsilon}{2}$-differential privacy.

Proof.

Computing the cost of a region that deviates from a uniform distribution, requires interactions with the histogram, and we must add sufficient noise to ensure differential privacy. We use the Laplace mechanism with $\epsilon_{1}$ to add random noise to the deviation cost. Note that the quadtree is just used as a search procedure to find a least-cost set of regions and we only reveal the leaves of quadtree (the uniform regions) for step 2 and the other regions are removed. Since there is no overlap among the leaves, we can use $\epsilon_{1}$ budget to compute the noisy value of uniformity cost for each region. Hence, the scale of noise is proportional to the sensitivity of regions, denoted by $\Delta cost$, instead of $O(4^{h+1})$ for the entire quadtree. To ensure privacy, we add a random $\nu_{i}\sim Lap(\Delta cost/\epsilon_{1})$ to the deviation cost of each region. Thus, the total cost of a region is defined in terms of the deviation cost and the error due to the perturbation noise:

[TABLE]

$\Delta cost$ shows the maximum change in a region’s cost when a trajectory is added to (or removed from) the data set. Assuming $k$ as the maximum length of a trajectory, in the worst case, a trajectory increases the count of $k$ cells in a region $p$, and we denote the $k$ cells by set $p^{\prime}$. We show $\Delta cost\leq 2k$.

Let $cost^{\prime}(p)$ be the region deviation after adding the trajectory. Assuming that the trajectory only affects the partition $p$ and $|p|\geq k$, we have:

[TABLE]

Line $(3)$ is derived using the triangle inequality and line $(4)$ follows from $k>0$. Using this result we have:

[TABLE]

Since a trajectory contributes to every level of the quadtree and thus would be counted multiple times in the overall trajectory count, we divide each trajectory contribution by the number of cells in which it is counted, i.e., its length. Using this approach, we reduce $k$ to $1$, i.e., its actual contribution to the overall trajectory count. This implies that $\Delta cost\leq 2$. We use the technique for computing the cost of regions.

In addition to the regions costs, in Algorithm 1 we compute the density of regions located in the leaves of the quadtree. Since computing the density requires accessing the original histogram $H$, we use the Laplace mechanism with the remaining privacy budget for step 1, $\epsilon_{2}$, to add $\nu^{\prime}\sim Lap(\frac{1}{\epsilon_{2}})$ noise to the density count for each region:

[TABLE]

The sensitivity of the density computation is $\leq 1$ as it follows the sensitivity of range queries (see Section II-B). Computing the region costs is $(\frac{\epsilon_{1}}{\Delta cost})$-differentially private and the density calculation is $(\epsilon_{2})$-differentially private. Given $\epsilon_{1}=\epsilon_{2}=\frac{\epsilon}{4}$ and using serial composition for differential privacy, Algorithm 1 is $(\frac{\epsilon}{4\Delta cost}+\frac{\epsilon}{4})\leq\frac{\epsilon}{2}$-differentially private. ∎

Accuracy. The accuracy is measured in terms of the difference between the resultant partitions and the true partitions for the given histogram.

Theorem III.2**.**

With probability at least $1-\alpha$, Algorithm 1 generates partitions with cost at most $TRU+\rho$, where $TRU$ is the cost of true partitions and $\rho=16nlog(n/\alpha)/\epsilon$.

Proof.

See Appendix A for the proof. ∎

Efficiency. Computing optimal partitions in a 2D space are computationally expensive. The complexity of a naive approach is $\Omega(n^{2})$ where $n$ is the number of cells in the histogram. Using the quadtree improves the computational complexity of partitioning to $O(nlog(n))$ as: (1) we partition the data space of a histogram with a fixed cell size and fixed counts, and (2) the structure is compatible with the trajectory data sets and each branch proceeds independently. While the achieved quadtree may not be balanced, the maximum height of the tree is of order $log(n)$ due to (1).

IV Step 2: Private Synthesizing Process

This section describes the second step of DQAM. The inputs are a set of queries $Q$, and the partitions $P=\{p_{1},...,p_{k}\}$ and densities $B=\{b_{1},...,b_{k}\}$ computed for the histogram $H$ in Section III. In this step, we aim to generate a synthetic histogram $H^{\prime}$ to be as close as possible to $H$. We believe this is the first method proposed that is both data- and query-aware for estimating an optimal spatial histogram $H^{\prime}$.

IV-A Synthesizing method

Our approach is based on the Multiplicative Weights Update Rule (MWUR) [7, 8], which is a query-aware method for estimating the distribution of histogram values. Starting from an initial uniform distribution, MWUR iteratively selects one query and updates the estimation by rescaling the values based on the query error. Intuitively, if we improve the estimated data for a query with the largest error, the error of other queries overlapping with this query will also be reduced. Hence, instead of considering all the given queries, MWUR only considers a small subset of queries to estimate the distribution of data with high utility for all queries. Based on this insight, we create a data- and query- aware update rule for privately synthesizing spatial histograms.

As MWUR is query-aware but not data-aware, we have to modify MWUR when synthesizing spatial histograms. Providing a data-aware approach incurs three challenges:

The cell counts in a region need to be updated/rescaled proportional to the density of that region. However, MWUR considers an equal weight/scale for all cells. Instead of assigning actual counts to each histogram cell, MWUR uses the relative proportion of counts (referred to as weight) in each cell. 2. 2.

MWUR only rescales the cell counts contributing to the selected queries. This rescaling strategy is problematic when regions partially overlap with the query area, and some of their cells are outside the query area. For such regions, MWUR only rescales the cells inside the query area and leaves the cells outside the query area unchanged. Thus, MWUR may degrade the utility of such regions especially those with uniform distribution when answering range queries. 3. 3.

The update function in MWUR cannot guarantee the local sequentiality of cells, which is a key property of spatial histograms. Ignoring this property can lead to an inconsistent output and in turn, reduces the utility.

In Section IV-B we address the first and second challenge by proposing a data- and query-aware update function. Given the uniform regions and their densities from step 1, our update function assigns a scaling factor to each region proportional to its density and preserves the uniformity of the updated regions. To address the third challenge, we need to ensure the consistency of the updated spatial histogram. A consistent spatial histogram guarantees that the answer of a range query monotonically increases with its query area when the larger range query includes smaller range queries. [6] developed an efficient approach to generate consistent output but the resulting histogram may not be optimal. In Section IV-C we discuss consistency for spatial histograms and develop a set of requirements to ensure consistency. We propose a highly efficient method to compute the optimal spatial histogram while ensuring consistency. The privacy budgets $\epsilon_{3}$ and $\epsilon_{4}$ are used for this step to ensure differential privacy. Algorithm 2 shows our data- and query-adaptive solution for estimating the true histogram.

IV-B Update function

Algorithm 2 shows, in each iteration $i\in T$, the Exponential mechanism consumes $\frac{T}{\epsilon_{3}}$ budget to privately select a query $q\in Q$ such that $q$ has worst error on the current estimation $H^{\prime}_{i}$. The error of selected query $q$ is then used to update the estimation, i.e.,

[TABLE]

$\nu_{i}\sim Lap(\frac{T}{\epsilon_{4}})$ is injected noise to the true answer of $q$.

In each iteration, MWUR reduces the estimation error by rescaling the weight of cell counts. Each cell count $H^{\prime}_{i}[j]$ intersecting the query area is multiplied by a scaling factor $s_{i}$ computed as:

[TABLE]

$q_{i}[j]$ is $j^{th}$ element value of query $q$ selected in iteration $i$. $q_{i}[j]=1$ if the cell with index $j$ is inside the query area and $q_{i}[j]=0$ otherwise. The scaling factor (i) only rescales the count of cells inside a query area and (ii) all the counts are rescaled equally. Such scaling function requires the original histogram $H$ to be uniformly distributed and fails if the densities in different regions have non-uniform distribution (see Section III-A).

We propose a scaling function that (i) updates all cells in a region if the region intersects the query area and (ii) rescales a cell count proportional to the density of region that the cell belongs to it. For a region $p\in P$ with density $b_{p}\in B$, our update function rescales the cell $p[j]$ with the following factor:

[TABLE]

where $q_{i}[j]=1$ if the region $p$ intersects $q$ and $q_{i}[j]=0$ otherwise. Note now the affected area by $q_{i}[j]$ is expanded beyond the query region. The query error $werr_{i}$ determines the magnitude and direction of update with respect to the selected query. $b_{p}$ amplifies the magnitude of update regarding the region density and $q_{i}[j]$ identifies the extent of update which is the union of query area and the regions intersected the query. According to the scaling factor, we define our update rule for iteration $i\in T$ as follows:

[TABLE]

where $n$ is the total number of trajectories in the histogram.

IV-C Optimal estimation

For a spatial histogram, it is key to preserve local sequentiality of cells, which means any change of a cell value cannot be performed independently without taking into account the adjacent cell values. As mentioned before, the edge values in spatial histograms capture the local sequentiality: an increase in a cell also leads to an increase in an adjacent cell if the edge value indicates that the represented trajectory (or more generally object) spans over both cells. Hence, local sequentiality imposes a constraint between each pair of adjacent faces $f_{i}$ and $f_{j}$ and their incident edge $e_{ij}$, which has to comply with the constraints $e_{ij}\leq f_{i}$ and $e_{ij}\leq f_{j}$ (a trip is always started in a cell per construction and might or might not cross an edge). This constraint is a key requirement for consistency in a spatial histogram as any change of data must satisfy it. Violating this constraint can even lead to a negative number of trajectories for a range query on a spatial histogram.

The update rule in Section IV-B does not consider this constraint when rescaling the face/edge values while computing a synthetic histogram which may lead to inconsistencies. Fig. 4 shows an example of the sequentiality violation where the query fully matches to the region. Since the update function rescales the values inside the query area it causes violations in edges placing on borders of the query. To guarantee the consistency of the synthesized histogram, we describe an approach using linear programming called consistent inference. Its aim is to compute the closest consistent histogram to the estimated histogram $H^{\prime}$ that satisfies the constraints.

Set of Consistency Constraints ($C$). Every edge count is at most equal to the minimum face count of its adjacent cells:

[TABLE]

Definition IV.1**.**

Consistent Inference. Let $T$ be the total number of iterations and $H^{\prime}$ be the estimated histogram in iteration $i\in T$. Given the set of constraints $\mathcal{C}$ between face and edge counts, the consistent inference function returns the histogram $H^{\prime\prime}=\{F^{\prime\prime},E^{\prime\prime}\}$ that satisfies the constraints in $\mathcal{C}$ while minimizing $\|H^{\prime\prime}-H^{\prime}\|$.

We formulate the consistent inference problem as follows:

[TABLE]

We use linear programming to solve this problem. Note that a constraint violation can only happen when the update rule scales down the counts or the query error in (5) is negative. Checking this condition in Algorithm 2 s3 significantly improves the efficiency of this algorithm.

IV-D Formal Guarantees

We evaluate Algorithm 2 in three key aspects: privacy, accuracy and efficiency.

Privacy. We show Algorithm 2 is $\frac{\epsilon}{2}$-differentially private.

Theorem IV.1**.**

Algorithm 2 satisfies $\frac{\epsilon}{2}$-differential privacy.

Proof.

Given $\epsilon_{3}=\epsilon_{4}=\frac{\epsilon}{4}$, the query selection and computing the error of selected query, each uses $\frac{\epsilon}{4T}$ of the privacy budget. Computing the optimal estimation is a post-processing and hence does not spend any privacy budget. The algorithm makes $T$ calls of the first two stages and according to the serial composition, the budget share accumulates in the calls which result $\frac{\epsilon}{2}$ privacy level for the algorithm. ∎

Accuracy. The accuracy is measured in terms of the maximum error bound in answering per $q\in Q$ using the estimated histogram $H^{\prime}$. The error is because of perturbation noise to ensure differential privacy and estimation error in rescaling the estimated counts in $H^{\prime}$ toward the true counts in $H$.

We show that the accuracy bound of Algorithm 2 meets the bounds with the state-of-the-art approaches in [7, 6]. However, this is the worst case bound. In our experiments, we show that defining the update rule based on regions and optimizing the estimation with respect to the consistency constraints significantly improves the output accuracy.

Theorem IV.2**.**

Given a spatial histogram $H$, for any $Q,T\in\mathbb{N}$, and $\epsilon>0$, with probability $\geq 1-2T/|Q|$, Algorithm 2 estimates $H^{\prime}$ such that

[TABLE]

Proof.

The proof is based on analyzing the effect of perturbation noise and the estimation error on the maximum error in answering queries. See appendix B for the proof. ∎

With the fact that $q(H)\in[0,n]$ for all $q\in Q$, we can reduce the bound to less than $n$ by choosing $T$ to be greater than $4log(|H|)$.

Efficiency. The computational cost of Algorithm 2 depends on the three steps in each iteration $i\in T$: query selection, updating the estimation, and computing the optimal estimation. For the query set $Q$, computing the score of each query $q\in Q$ has $O(|H|)$ cost where $|H|=r\times c$ is the total number of cells in $H$. Thus, the total cost of query selection is $O(|H||Q|)$. Updating the estimation takes $O(|H|)$ as the cell counts in regions overlapping with the query should be rescaled by the scaling factor. The linear programming solution in the last step has the worst case cost of $O(|H|^{3.5})$ [11]. The three steps are called for $T$ times that results in the complexity of $O(T(|H||Q|+|H|^{3.5}))$ for the algorithm.

In practice, the total number of queries selected through the iterations (one per iteration $i\in T$), call $Q^{\prime}$, is much less than the $|Q|$, i.e., $|Q^{\prime}|\ll|Q|$. The $T$ is also chosen to be a small constant such that $T\cong log(|H|)$ due to the accuracy bound above. These two small constants have a negligible effect on the first part of the complexity cost. In addition, using the violation check after an update can significantly improve the complexity of the algorithm as we do not need to call linear programming function.

V Experimental Evaluation

This section evaluates the performance of DQAM in different settings explored by recent works including privacy budgets, query sets, data sets, the entire world size (i.e., (1) number of records in the trajectory data set, and (2) number of cells in the histogram), and running time. We contrast the quality of generated histogram with the output of recently proposed algorithms in terms of accuracy and utility.

V-A Experimental setup

Considering the number of cells (rows times columns) in a histogram as its resolution, we set the resolution of a histogram as a power of $2$. Depending on the speed limit of moving objects, each data set may need a different level of resolution for the histogram. In a data set with $50km/h$ limit and trajectory sampling rate of 3 GPS coordinates/sec, a histogram with resolution $8$ aggregates the trajectory counts into $2^{8}\times 2^{8}$ cells with $39m^{2}$ per cell which is fairly accurate for capturing movements in metropolitan areas. Higher resolutions can more accurately capture the movements but are computationally more expensive. Unless otherwise specified, we map the data sets into histograms of resolution $8$. We choose the partitioning threshold $\delta=4/\epsilon_{1}^{2}$ in Algorithm 1 to cover the variance of total noise added to the costs of a partition and its potential children. This choice of $\delta$ prioritizes high uniformity in partitioning which leads to highly uniform partitions. We found the computed results stable with this setting for $\delta$. The number of iterations $T$ is chosen from the set $\{10,12,14,16,18,20\}$. For each chosen $T$ value, we run the algorithm $5$ times and compute the mean error. Finally, we report the lowest error.

Data sets. One real and one synthetic data set are used in our experiments. The real data set is Porto-Taxi containing the trajectory of taxi-cabs in the city of Porto, Portugal [19]. It has about $1.7$ million trajectories of various lengths from $443$ taxi-cabs in one year. Porto-Taxi is a large data set with a highly skewed distribution of trajectories going toward the city center that mapping it to a spatial histogram, results in a histogram with many sparse areas. It is important for the algorithm to be able to capture the distribution information and utilize it by privately estimating the true distribution. The synthetic data set is called Melb-Car which is generated by Minnesota Web-based Traffic Generator (MNTG) [17]. We simulated $1,000,000$ vehicle trajectories for $20$ time units in the city of Melbourne, Australia. Melb-Car represents a data set with a uniform distribution where the counts are almost equal in all histogram cells. The two data sets are used to measure the effect of data set distribution on the quality of estimation. To show the effect of data set size $|D|$ on the estimation quality, we use three samples of each data set in the experiments, $|D|\in\{1000,$100,000$,$1,000,000$\}$. We used a square of $100$km2 as the spatial geometry to map the data set into a histogram. Unless otherwise specified, a data set of $1000$ trajectories from Porto-Taxi is used for the experiment. We note such a small data set with non-uniform distribution results in a sparse histogram and provides the most difficult setting for the algorithm to estimate the true distribution of data. Our results are of high quality even in this hard setting.

Queries. We use queries in the experiments that are uniformly sampled from the domain of histogram $H$. Given $X=[x,x^{\prime}]$ and $Y=[y,y^{\prime}]$ as domains of columns and rows in $H$ respectively, a query is generated as a rectangle $([x_{1},x_{2}],[y_{1},y_{2}])$ where $x_{1},x_{2}\in X$ and $x_{1}\leq x_{2}$; $y_{1},y_{2}\in Y$ and $y_{1}\leq y_{2}$. Unless otherwise specified, the default size $|Q|=$16,000$$ is used for the query set.

Measures. We evaluate the accuracy of estimation in terms of the error in answering queries and measure it using the average $L_{1}$ error per query. In addition, we expect the estimated histogram have high utility in answering arbitrary queries from the class of range queries. The utility of estimated histogram $H^{\prime}$ is evaluated as the difference between its distribution $H^{\prime}/n$ and the distribution of original histogram $H/n$:

[TABLE]

where $n$ is the total number of trajectories in the given data set. This measure is also known as Kullback-Leibler Divergence (KLD) and is used in prior work [7, 6] to evaluate the utility of estimation. This enables us to compare our approach to previous work.

Baselines. We compare the performance of DQAM with four existing works: Laplace mechanism (LM) [5] as a naive approach, MWEM [7], DAWA [13] and PriSH [6]. Table II contrasts the baselines and DQAM in terms of technical features, i.e., (1) utilizing data properties (data-aware), (2) utilizing query properties (query-aware), and and the output features, i.e., (3) preserving consistency or local sequentiality (greedily consistent) as well as (4) Optimizing the output while ensuring consistency (optimally consistent). As the table shows, LM and MWEM only consider the query properties to identify the scale of noise. LM is a naive approach in which the noise is directly scaled based on the dependency of queries. In trajectory data, the sensitivity is equivalent to the number of queries which is considerably large. MWEM, however, selects a few numbers of queries which significantly reduces the scale of added noise. DAWA, is data-aware and query-aware similar to DQAM but does not consider the local sequentiality in trajectories (see Section VI for more details). Hence, DAWA fails to preserve consistency in the output. On the other hand, PriSH maintains the consistency of data. However, it only takes advantage of query properties in scaling the noise and generating the output histogram. Additionally, Prish fails to generate the optimal histogram with respect to the consistency requirements (see Section VI for more details). DQAM is the only mechanism that is data- and query-aware, and generates a histogram that is optimal with respect to consistency preservation.

Since MWEM and DAWA are designed for 1D histograms with no local sequentiality, we only use the face component to evaluate the algorithms. For these algorithms, we convert the 2D histogram to 1D using Hilbert curves. The number of iterations in PriSH and MWEM are chosen as explained for DQAM, and $\epsilon$ is split evenly between the query selection and the update function in the algorithms. According to [13], we consider the ratio $r=0.5$ for dividing privacy budget $\epsilon$ between the data-aware and query-aware steps in DAWA. In our evaluations, we do not report results of DQAM and PriSH algorithms without consistency. Note that LM, MWEM and DAWA do not guarantee the consistency of generated histogram. Since we only consider the face component to evaluate these mechanisms no penalty is applied for the inconsistency of output. The corresponding results of LM, MWEM and DAWA are depicted by dashed lines in graphs to indicate this limitation.

V-B Accuracy Evaluations

Data size and distribution. Fig. 5 evaluates the accuracy of DQAM versus existing works in different data set sizes where results are depicted in log scale. For larger datasets, the counts in the corresponding histogram are larger and hence, less sensitive to noise. That is the algorithms usually have higher accuracy for larger data sets. As the results show, DQAM achieves mostly the lowest error specifically for smaller $\epsilon$s. This advantage is more significant on Porto-Taxi data set where DQAM accurately captures the sparse/dense regions and estimates the counts properly. In contrast, PriSH has an unstable accuracy across different $\epsilon$ values as the heuristic strategy in PriSH fails to efficiently estimate histogram in the presence of noise.MWEM and DAWA work good on uniform distribution in Melb-Car, but their accuracy considerably decreases on non-uniform distribution. The generated histogram, though, is not valid because of inconsistency in data. The error of generated histogram by LM is always considerably higher than other mechanisms due to naively scaling the added noise. Table III reports the standard deviation of accuracy for algorithms for different $\epsilon$ values. As expected, the standard deviation of LM is always high, and this naive approach does not show a stable behaviour. DQAM has the least standard deviation. In contrast with PriSH, DQAM employs the optimal consistent estimation resulting in smaller and more stable error which decreases for larger $\epsilon$s. Inconsistency of generated histograms by MWEM and DAWA result in large variations in the error.

Size of query set. As the set of given queries are randomly generated, the queries may have different shapes. This diversity helps the query selection in DQAM to capture the information from the spatial histogram efficiently. We evaluate the effect of changing the query set size on the error of estimated histogram. Fig. 6 depicts the error of DQAM versus baseline methods. DQAM often achieves higher accuracy comparing to the baselines. Furthermore, the error is almost stable and slightly increases in larger query sets. This behaviour agrees with the Theorem IV.2 in that the size of the query set has a logarithmic effect on the worst case error.

Resolution of histogram. Fig. 7 investigates the effect of changing the histogram resolution on the quality of estimation in DQAM. As it was expected, higher resolutions cause higher errors. However, DQAM works well in different resolutions comparing to the other algorithms. The significant increment in the error of LM, MWEM and DAWA is due to accumulation of error caused by local sequentiality violations.

As mentioned above, the reported results are the average value of $5$ independent runs of an algorithm. The ratio of accuracy improvement achieved by DQAM comparing to the competing algorithms changes in different settings. The highest ratio achieved in our experiments was $7.4$ for $|D|=1000$, $|Q|=8000$, resolution $8$ and $\epsilon=0.01$ comparing the average error of DQAM versus PriSH shown in Fig. 6: PriSH/DQAM $=14.0654/1.901=7.3989$. While the ratio of DQAM accuracy versus LM, MWEM and DAWA is significantly higher (minimum $\geq 60$ for this setting), we did not report it as the output of these algorithms are not consistent due to local sequentiality violations.

V-C Utility Evaluation

The utility of estimated histogram shows its quality in representing the original histogram. Utility is in fact a measure of quality of generated histogram in answering arbitrary range queries, the ones not included in the query set or not selected in the mechanism process. As mentioned before, We measure the utility of estimated histograms by $KLD$ which is used in prior works. We use this measure to show employing both data and query information in DQAM has significant effect in improving the algorithm in detecting sparse/dense regions. Fig. 8 depicts that DQAM perfectly captures the data distribution properties in both uniform and non-uniform distributions. For clarity in presentation, we focus on DQAM, PriSH and MWEM in this experiment.

V-D Efficiency Evaluation

Fig. 9 evaluates the running time of DQAM on small to large samples of Porto-Taxi and Melb-Car datasets with $\epsilon=0.1$ and a fixed resolution of histogram. The partitioning (step 1 in DQAM) takes less than a second in all settings, and the synthesizing process (step 2 in DQAM) mostly contributes to the running time. The total running time for large samples of Porto-Taxi and Melb-Car are $\approx 10$ and $\approx 8$ seconds, respectively. In small samples where the histogram is very sparse, some updates in the synthesizing process cause consistency violations which, in turn, cause the second step to take more time to compute an optimal histogram satisfying consistency. Note that a violation should be corrected as soon as it happened in the synthesizing process as it can affect the updates in next iterations. Indeed, the accumulation of consistency violations leads the solution to deviate in a wrong direction that it may be far from the optimal point that causes the mechanism to take a long time to compute the optimal estimation in the feasible space of the problem. The total running time for small datasets is $\approx 60$ minutes for Porto-Taxi and $\approx 25$ minutes for Melb-Car. According to the results our mechanism is eminently practical and scalable.

VI Related Works

Studies on human movement patterns [4, 25] show the trajectories are highly identifiable and that coarsening or aggregation barely reduce their uniqueness. Montjoye et al. [4] showed that even after coarsening trajectories, $95$% of individuals could be uniquely identified with only $4$ location points from their trajectories. Recently, Xu et al. [25] studied the effectiveness of aggregation in anonymizing trajectories in large-scale data sets and $91$% of trajectories could be uniquely reconstructed. Hence, private mechanisms with a strong privacy guarantee are needed to publish trajectories.

To address this problem, data-aware mechanisms introduce $\epsilon$-differential privacy to the aggregation techniques [2, 10] and reconstruct the trajectory data set from the aggregated data. Even though the mechanisms capture the local sequentiality in trajectories, they are not ideal for efficiently answering range queries. Chen et al. [2] use a Markov process to model local sequentiality and construct a hierarchical structure such as a prefix tree by grouping the trajectories with the same prefix. Subsequent works [1, 10] use most frequent substrings to get higher counts in the leaves and to achieve a better utility. However, a recent work [20] showed that the mechanisms fail to preserve trajectory properties such as their length, and thus cannot guarantee the accurate distribution of trajectories. Even the use of a hierarchical structure, i.e., the constructed noisy prefix tree, is computationally expensive for range queries. Given a prefix tree of height $l>0$ with a branching factor $a>1$, answering a range query has a computational cost of $O(a^{l+1})$. Our mechanism, however, publishes a 2D histogram that maintains the local sequentiality and costs $O(k)$ on range queries where $k$ is the size of query area.

Monreale et al. [18] define the dependency between cells instead of points by mapping the trajectories to a grid to count the movement frequencies between the adjacent cells. However, a frequency vector only maintains the number of transitions for a group of observations without information about the spatial adjacency of two vectors. This information is crucial for a range query as it counts the vectors overlapping the query area. In our work, we use a spatial histogram that is a grid but captures both the cell counts and the spatial adjacency of trajectories.

Related work in coarsening the data space are in two categories: data-aware mechanisms [3, 21] and query-aware mechanisms [5, 7, 13]. The data-aware mechanisms use some data properties to adaptively partition the data space and optimize the error of published data in answering range queries. However, the models are developed for location points rather than trajectories. Cormode et al. [3] build a hybrid tree in which a kd-tree partitions the space into coarse density regions and a quadtree maps the dense areas into finer resolutions. Qardaji et al. [21], generate partitions to get uniform regions of data. The partitioning condition is based on the density of location points in regions. The data-aware partitioning techniques make use of information about the data distribution to enhance accuracy. However, the mechanisms are not able to capture local sequentiality of location points which is crucial in representing trajectories. The query-aware mechanisms, on the other hand, take a set of range queries to get information about the data distribution. Laplace mechanism [5] is a naive query-aware model that adds an independent noise to the answer of each query in the query set. Hardt et al. [7] developed MWEM that optimizes the accuracy of published data for the class of range queries. By laying a grid with a fixed cell size on the data space, MWEM iteratively learns the density of cells using the noisy answer of queries. Li et al. [13] proposed a data- and query-aware (DAWA) which employs data properties as well as query properties in publishing data under differential privacy. A comparative analysis of space decomposition mechanisms by Hay et al. [9] showed MWEM and DAWA achieve high accuracy in answering range queries on location data, not trajectories. Since the mechanisms are not able to capture the local dependency between locations which is crucial in studying trajectories. In this paper, our proposed mechanism uses spatial histograms in which the edge component maintains the local sequentiality and the learning process ensures preserving this property in the published spatial histogram.

Recently, Ghane et al. [6] developed a mechanism named Private Spatial Histogram PriSH for range queries on trajectories. PriSH publishes a synthetic spatial histogram under $\epsilon$-differential privacy. It is a query-aware mechanism that extends the idea of MWEM in [7]. PriSH, takes a spatial histogram and a query set as input and utilizes the correlation between the queries to estimate the distribution of the original histogram privately. To maintain the consistency in the histogram, PriSH uses a heuristic approach that may result in a histogram far from the original spatial histogram. The reason lies in the heuristic approach that locally ensures consistency and may lead to overcorrection. In this paper, we propose a data- and query- aware mechanism that utilizes the trajectories density in different regions as well as the given queries correlation to estimate the optimal spatial histogram with significantly higher utility. Our query-aware strategy employs a linear programming approach to provide a guarantee on optimally consistent histogram which leads to significant improvement in the utility of results.

VII Conclusion

We proposed DQAM, a data and query-aware mechanism for privately publishing spatial histograms and answering range queries on trajectories. The first step of DQAM identifies disjoint uniform regions in a histogram and computes the density of each histogram region. It uses a quadtree as it enables efficient space partitioning. Based on the quadtree, we can measure the uniformity for each step of the partitioning. The partitions and densities are then used in the second step of DQAM for learning the true distribution of histogram. The learning process uses a given query set as well as the partitions and densities from the first step to estimate the true distribution and to improve the quality of estimation in each iteration. Our experiments show that DQAM significantly improves the accuracy over the state of the art methods, specifically on data sets with a non-uniform distribution while ensuring consistency, i.e., local sequentiality. Our experimental results show that we have improved utility by a factor 7.4 (see Fig. 6). In the future, we will focus on differentially private trajectory representations that retain high data utility and thus can be used for any spatial query type as current techniques cannot achieve the required data utility.

Appendix A Proof of Theorem III.2

Let $TRU$ be the cost of true partitions and $L_{\rho}$ be the set of all possible partitions that for $\textbf{P}\in L_{\rho}$ the total cost of partitioning is $cost(\textbf{P})\leq TRU+\rho$. $\tilde{L_{\alpha}}$ represents the complement of set $L_{\alpha}$. Let $\nu$ be the Laplace random noise with scale $b=\Delta cost/\epsilon_{1}$ added to each region $p\in\textbf{P}$. Given $|H|=n$ as the number of cells in the histogram, if $\mathinner{\!\left\lVert\nu\right\rVert}<\frac{\rho}{2n}$, the total noise added to the partition cost is $n.\frac{\rho}{2n}=\frac{\rho}{2}$ because a partitioning can have at most $n$ regions. This means with such noise still no partitioning in $\tilde{L_{\alpha}}$ will be generated. Given Algorithm 1 as $\mathcal{A}$, we have:

[TABLE]

The line $(7)$ is achieved by union bound. The line $(8)$ follows the definition of Exponential mechanism with scale $b$ because the absolute value of a random variable from Laplace distribution has an exponential distribution. Now, we can compute the $\rho\geq 2.b.n.log(n/\alpha)=\frac{4.n.log(n/\alpha)}{\epsilon_{1}}=\frac{16.n.log(n/\alpha)}{\epsilon}$

Appendix B Proof of Theorem IV.2

The proof follows the analyze of accuracy in [7]. We show that our new update rule and later the linear programming do not affect the worst case bound of accuracy. Without further amendments, the proof in [7] carries over from 1D to 2D for the spatial histogram in our work. It first computes a bound on the perturbation noise regarding the query error. Second, it uses the relative entropy to bounds the estimation error imposed by the update rule. Last, it combines the two bounds to achieve the worst case error bound on the accuracy. Our improvements affect the second step of analysis and we show it does not change the bound of estimation error.

Our update rule rescales the count of each cell $x\in H^{\prime}$ proportional to the density of region that $x$ belongs to that, i.e., $b_{x}\in B$. After scaling the counts in each round, the optimal estimator multiplies the cell counts by a multiplier $\alpha_{x}$ where $\alpha_{x}\geq 0$. Considering a maximum multiplier in each round, we define $\alpha$ as the maximum bound on the multipliers across all rounds, i.e.,

[TABLE]

According to the definition of relative entropy, the difference between the $H$ and $H^{\prime}$ after each update round improves as:

[TABLE]

where the ratio $c(H^{\prime}_{i}[x])/c(H^{\prime}_{i-1}[x])$ can be written as $\frac{\alpha}{\beta_{i}}sp_{i}$ such that $sp_{i}=exp(q_{i}[x]\eta_{i}b_{x})$ with $\eta_{i}=(m_{i}-q_{i}(H^{\prime}_{i-1})/2n$ and $b_{x}\in B$ is the density of region that the cell $x$ belongs to it. $\beta_{i}$ is the re-normalization factor. As the cells may have different density coefficients, we consider $\textbf{b}=\{\textbf{b}\in B\mid\forall b_{i}\in B,\textbf{b}\geq b_{i}\}$. Expanding the equation gives:

[TABLE]

Using the Taylor expansion of $sp_{i}$ and $q_{i}[x]^{2}\leq 1$, the $\beta_{i}$ bound will be:

[TABLE]

Inserting this bound to the equality 11 and that $log(1+x)\leq x$ and $\textbf{b}\leq 1$, cancels $log(\alpha)$ and gives:

[TABLE]

By introducing the definition of $\eta_{i}$ to this result, it gives:

[TABLE]

which is the same bound on the estimation error in [7]. The remaining of the proof follows the accuracy analysis in [7].