Proceedings on Privacy Enhancing Technologies ..; .. (..):1–22

Takao Murakami*, Koki Hamada, Yusuke Kawamoto, and Takuma Hatano

Privacy-Preserving Multiple Tensor
Factorization for Synthesizing Large-Scale
Location Traces with Cluster-Specific Features
Abstract: With the widespread use of LBSs (Locationbased Services), synthesizing location traces plays an increasingly important role in analyzing spatial big data
while protecting user privacy. In particular, a synthetic
trace that preserves a feature specific to a cluster ofusers (e.g., those who commute by train, those who
go shopping) is important for various geo-data analysis
tasks and for providing a synthetic location dataset. Although location synthesizers have been widely studied,existing synthesizers do not provide sufficient utility,
privacy, or scalability, hence are not practical for largescale location traces. To overcome this issue, we propose
a novel location synthesizer called PPMTF (PrivacyPreserving Multiple Tensor Factorization). We model
various statistical features of the original traces by atransition-count tensor and a visit-count tensor. We factorize these two tensors simultaneously via multiple tensor factorization, and train factor matrices via posterior
sampling. Then we synthesize traces from reconstructedtensors, and perform a plausible deniability test for asynthetic trace. We comprehensively evaluate PPMTF
using two datasets. Our experimental results show that
PPMTF preserves various statistical features including cluster-specific features, protects user privacy, andsynthesizes large-scale location traces in practical time.
PPMTF also significantly outperforms the state-of-theart methods in terms of utility and scalability at thesame level of privacy.
Keywords: location privacy, location synthesizer,
cluster-specific feature, multiple tensor factorization
DOI Editor to enter DOI
Received ..; revised ..; accepted ...

*Corresponding Author: Takao Murakami: AIST, Email: takao-murakami at aist.go.jp
Koki Hamada: NTT/RIKEN E-mail: hamada.koki atlab.ntt.co.jp
Yusuke Kawamoto: AIST, E-mail: yusuke.kawamoto ataist.go.jp
Takuma Hatano: NSSOL, E-mail: hatano.takuma.hq2 at
jp.nssol.nipponsteel.com

1 Introduction
LBSs (Location-based Services) have been used in avariety of applications such as POI (Point-of-Interest)
search, route finding, and geo-social networking. Consequently, numerous location traces (time-series location trails) have been collected into the LBS provider.
The LBS provider can provide these location traces(also called spatial big data [60]) to a third party (ordata analyst) to perform various geo-data analysis tasks;e.g., finding popular POIs [75], semantic annotationof POIs [19, 70], modeling human mobility patterns[17, 40, 42, 63], and road map inference [5, 41].
Although such geo-data analysis is important for industry and society, some important privacy issues arise.
For example, users’ sensitive locations (e.g., homes, hospitals), profiles (e.g., age, profession) [33, 44, 73], activities (e.g., sleeping, shopping) [39, 73], and social relationships [6, 24] can be estimated from traces.
Synthesizing location traces [8, 15, 32, 36, 64, 72]
is one of the most promising approaches to performgeo-data analysis while protecting user privacy. This approach first trains a generative model from the original
traces (referred to as training traces). Then it generates
synthetic traces (or fake traces) using the trained generative model. The synthetic traces preserve some statistical features (e.g., population distribution, transition
matrix) of the original traces because these features aremodeled by the generative model. Consequently, basedon the synthetic traces, a data analyst can perform thevarious geo-data analysis tasks explained above.
In particular, a synthetic trace that preserves a feature specific to a cluster of users who exhibit similar
behaviors (e.g., those who commute by car, those who
often go to malls) is important for tasks such as semantic annotation of POIs [19, 70], modeling human mobility patterns [17, 40, 42, 63], and road map inference[5, 41]. The cluster-specific features are also necessary
for providing a synthetic dataset for research [30, 52] oranonymization competitions [2]. In addition to preserving various statistical features, the synthetic traces are(ideally) designed to protect privacy of users who pro-Privacy-Preserving Multiple Tensor Factorization for Synthesizing Large-Scale Location Traces with Cluster-Specific Features

vide the original traces from a possibly malicious data
analyst or any others who obtain the synthetic traces.
Ideally, a location synthesizer should satisfy the following three features: (i) high utility: it synthesizes
traces that preserve various statistical features of theoriginal traces; (ii) high privacy: it protect privacy ofusers who provide the original traces; (iii) high scalability: it generates numerous traces within an acceptable time; e.g., within days or weeks at most. All ofthese features are necessary for spatial big data analysis or providing a large-scale synthetic dataset.
Although many location synthesizers [8, 12, 13, 15,
28, 32, 36, 64, 72] have been studied, none of them aresatisfactory in terms of all three features:
Related Work. Location privacy has been widely
studied ([11, 27, 37, 56] presents related surveys) andsynthesizing location traces is promising in terms of
geo-data analysis and providing a dataset, as explained
above. Although location synthesizers have been widely
studied for over a decade, Bindschaedler and Shokri [8]showed that most of them (e.g., [15, 32, 36, 64, 72])
do not satisfactorily preserve statistical features (especially, semantic features of human mobility, e.g., “manypeople spend night at home”), and do not provide highutility.
A synthetic location traces generator in [8] (denotedby SGLT) is a state-of-the-art location synthesizer. SGLT
first trains semantic clusters by grouping semanticallysimilar locations (e.g., homes, offices, and malls) basedon training traces. Then it generates a synthetic tracefrom a training trace by replacing each location with all
locations in the same cluster and then sampling a trace
via the Viterbi algorithm. Bindschaedler and Shokri [8]showed that SGLT preserves semantic features explained
above and therefore provides high utility.
However, SGLT presents issues of scalability, whichis crucially important for spatial big data analysis.
Specifically, the running time of semantic clustering in
SGLT is quadratic in the number of training users andcubic in the number of locations. Consequently, SGLT
cannot be used for generating large-scale traces. For example, we show that when the numbers of users and locations are about 200000 and 1000, respectively, SGLT
would require over four years to execute even by using
1000 nodes of a supercomputer in parallel.
Bindschaedler et al. [9] proposed a synthetic data
generator (denoted by SGD) for any kind of data using
a dependency graph. However, SGD was not applied tolocation traces, and its effectiveness for traces was unclear. We apply SGD to location traces, and show that it

2

cannot preserve cluster-specific features (hence cannot
provide high utility) while keeping high privacy. Similarly, the location synthesizers in [12, 13, 28] generate
traces only based on parameters common to all users,
and hence do not preserve cluster-specific features.
Our Contributions. In this paper, we propose a novel
location synthesizer called PPMTF (Privacy-Preserving
Multiple Tensor Factorization), which has high utility,
privacy, and scalability. Our contributions are as follows:– We propose PPMTF for synthesizing traces.
PPMTF models statistical features of training
traces, including cluster-specific features, by two
tensors: a transition-count tensor and visit-counttensor. The transition-count tensor includes a transition matrix for each user, and the visit-count tensor includes a time-dependent histogram of visited
locations for each user. PPMTF simultaneously factorizes the two tensors via MTF (Multiple Tensor
Factorization) [35, 65], and trains factor matrices(parameters in our generative model) via posterior
sampling [67]. Then it synthesizes traces from reconstructed tensors, and performs the PD (Plausible
Deniability) test [9] to protect user privacy. Technically, this work is the first to propose MTF in a
privacy preserving way, to our knowledge.
– We comprehensively show that the proposed
method (denoted by PPMTF) provides high utility,
privacy, and scalability (for details, see below).
Regarding utility, we show that PPMTF preserves all of
the following statistical features.
(a) Time-Dependent Population Distribution. The
population distribution (i.e., distribution of visited locations) is a key feature to find popular POIs [75]. It
can also be used to provide information about the number of visitors at a specific POI [29]. The population
distribution is inherently time-dependent. For example,
restaurants have two peak times corresponding to lunch
and dinner periods [70].
(b) Transition Matrix. The transition matrix is amain feature for modeling human movement patterns[42, 63]. It is used for predicting the next POI [63] orrecommending POIs [42].
(c) Distribution of Visit-Fractions. A distributionof visit-fractions (or visit-counts) is a key feature forsemantic annotation of POIs [19, 70]. For example, [19]reports that many people spend 60% of the time at theirhome and 20% of the time at work/school. [70] reportsthat most users visit a hotel only once, whereas 5% ofusers visit a restaurant more than ten times.

Privacy-Preserving Multiple Tensor Factorization for Synthesizing Large-Scale Location Traces with Cluster-Specific Features(d) Cluster-Specific Population Distribution. At
an individual level, a location distribution differs fromuser to user, and forms some clusters; e.g., those wholive in Manhattan, those who commute by car, and those
who often visit malls. The population distribution forsuch a cluster is useful for modeling human location
patterns [17, 40], road map inference [5, 41], and smart
cities [17].
We show that SGD does not consider cluster-specificfeatures in a practical setting (similarly, [12, 13, 28]
do not preserve cluster-specific features), and therefore
provides neither (c) nor (d). In contrast, we show that
PPMTF provides all of (a)-(d). Moreover, PPMTF automatically finds user clusters in (d); i.e., manual clustering is not necessary. Note that user clustering is very
challenging because it must be done in a privacy preserving manner (otherwise, user clusters may reveal information about users who provide the original traces).
Regarding privacy, there are two possible scenarios
about parameters of the generative model: (i) the parameters are made public and (ii) the parameters are
kept secret (or discarded after synthesizing traces) andonly synthetic traces are made public. We assume scenario (ii) in the same way as [8]. In this scenario, PPMTF
provides PD (Plausible Deniability) in [9] for a synthetic
trace. Here we use PD because both SGLT [8] and SGD[9] use PD as a privacy metric (and others [12, 13, 28] do
not preserve cluster-specific features). In other words,
we can evaluate how much PPMTF advances thestate-of-the-art in terms of utility and scalability at the same level of privacy. We also empiricallyshow that PPMTF can prevent re-identification (or deanonymization) attacks [26, 47, 62] and membership inference attacks [31, 61] in scenario (ii). One limitationis that PPMTF does not guarantee privacy in scenario(i). We clarify this issue at the end of Section 1.
Regarding scalability, for a larger number |U| oftraining users and a larger number |X | of locations,
PPMTF’s time complexity O(|U||X |2 ) is much smaller
than SGLT’s complexity O(|U|2 |X |3 ). Bindschaedler andShokri [8] evaluated SGLT using training traces of only
30 users. In this paper, we use the Foursquare dataset in[68] (we use six cities; 448839 training users in total) andshow that PPMTF generates the corresponding traces
within 60 hours (about 106 times faster than SGLT) by
using one node of a supercomputer. PPMTF can also
deal with traces of a million users.
In summary, PPMTF is the first to provide all of
the utility in terms of (a)-(d), privacy, and scalability
to our knowledge. We implemented PPMTF with C++,3

and published it as open-source software [1]. PPMTF
was also used as a part of the location synthesizer toprovide a dataset for an anonymization competition [2].
Limitations. Our results would be stronger if user privacy was protected even when we published the parameters of the generative model; i.e., scenario (i). However,
PPMTF does not guarantee meaningful privacy in this
scenario. Specifically, in Appendix F, we use DP (Differential Privacy) [21, 22] as a privacy metric in scenario(i), and show that the privacy budget ε in DP needs tobe very large to achieve high utility. For example, if we
consider neighboring data sets that differ in one trace,then ε needs to be larger than 2×104 (which guarantees
no meaningful privacy) to achieve high utility. Even if
we consider neighboring data sets that differ in a single location (rather than one trace), ε = 45.6 or more.
We also explain the reason that a small ε is difficult inAppendix F. We leave providing strong privacy guarantees in scenario (i) as future work. In Section 5, we also
discuss future research directions towards this scenario.

2 Preliminaries
2.1 Notations
Let N, Z≥0 , R, and R≥0 be the set of natural numbers, non-negative integers, real numbers, and nonnegative real numbers, respectively. For n ∈ N, let[n] = {1, 2, · · · , n}. For a finite set Z, let Z ∗ be theset of all finite sequences of elements of Z. Let P(Z) be
the power set of Z.
We discretize locations by dividing the whole map
into distinct regions or by extracting POIs. Let X be a
finite set of discretized locations (i.e., regions or POIs).
Let xi ∈ X be the i-th location. We also discretize time
into time instants (e.g., by rounding down minutes toa multiple of 20, as in Figure 1), and represent a time
instant as a natural number. Let T ⊂ N be a finite setof time instants under consideration.
In addition to the time instant, we introduce a timeslot as a time resolution in geo-data analysis; e.g., if wewant to compute the time-dependent population distribution for every hour, then the length of each time slot
is one hour. We represent a time slot as a set of time
instants. Formally, let L ⊆ P(T ) be a finite set of timeslots, and li ∈ L be the i-th time slot. Figure 1 shows anexample of time slots, where l1 = {1, 2, 3}, l2 = {4, 5, 6},l3 = {7, 8, 9}, and L = {l1 , l2 , l3 }. The time slot cancomprise either one time instant or multiple time in-Privacy-Preserving Multiple Tensor Factorization for Synthesizing Large-Scale Location Traces with Cluster-Specific Features

User
u1

s1

x2

x3

u2

s2

x1

u3

s3

x5

x4

Hour (AM) 7
Time Instant 1
Time Slot

2

x4
x1

l1

3

Training Tracesx3 x4 x5
x2

x4

x3
x5

x1

x3

x5

x3

8
4

5

6

9
7

l2

x5

8

x3

9

4

of un , whereas it preserves a lot of features specific toclusters un belongs to. Therefore, there is a trade-off between the cluster-specific features and user privacy. In
Appendix C, we show an example of Mn in SGD [9].
In Appendix A, we also show tables summarizing
the basic notations and abbreviations.

l3Fig. 1. Training traces (|U | = 3, |X | = 5, |T | = 9, |L| = 3).
Missing events are marked with gray.

stants (as in Figure 1). The time slot can also compriseseparated time instants; e.g., if we set the interval between two time instants to 1 hour, and want to average
the population distribution for every two hours over twodays, then l1 = {1, 2, 25, 26}, l2 = {3, 4, 27, 28}, · · · , l12 ={23, 24, 47, 48}, and L = {l1 , · · · l12 }.
Next we formally define traces as described below.
We refer to a pair of a location and a time instant as an
event, and denote the set of all events by E = X × T .
Let U be a finite set of all training users, and un ∈ U
be the n-th training user. Then we define each trace
as a pair of a user and a finite sequence of events, anddenote the set of all traces by R = U × E ∗ . Each trace
may be missing some events. Without loss of generality,
we assume that each training user has provided a single
training trace (if a user provides multiple temporallyseparated traces, we can concatenate them into a single
trace by regarding events between the traces as missing).
Let S ⊆ R be the finite set of all training traces, andsn ∈ S be the n-th training trace (i.e., training trace ofun ). In Figure 1, s1 = (u1 , (x2 , 1), (x3 , 2), (x4 , 3), (x3 , 5),(x4 , 6), (x5 , 7), (x5 , 9)) and S = {s1 , s2 , s3 }.
We train parameters of a generative model (e.g., semantic clusters in SGLT [8], factor matrices in PPMTF)
from training traces, and use the model to synthesizea trace. Since we want to preserve cluster-specific features, we assume a type of generative model in [8, 9]
as described below. Let y ∈ R be a synthetic trace.
For n ∈ [|U|], let Mn be a generative model of user unthat outputs a synthetic trace y ∈ R with probabilityp(y = Mn ). Mn is designed so that the synthetic trace y(somewhat) resembles the training trace sn of un , while
protecting the privacy of un . Let M be a probabilistic generative model that, given a user index n ∈ [|U|]as input, outputs a synthetic trace y ∈ R producedby Mn ; i.e., p(y = M(n)) = p(y = Mn ). M consistsof M1 , · · · , M|U | , and the parameters of M1 , · · · , M|U |are trained from training traces S. A synthetic tracey that resembles sn too much can violate the privacy

2.2 Privacy Metric
We explain PD (Plausible Deniability) [8, 9] as a privacy metric. The notion of PD was originally introduced
by Bindschaedler and Shokri [8] to quantify how well a
trace y synthesized from a generative model M provides
privacy for an input user un . However, PD in [8] was defined using a semantic distance between traces, and its
relation with DP was unclear. Later, Bindschaedler et
al. [9] modified PD to clarify the relation between PD
and DP. In this paper, we use PD in [9]:Definition 1 ((k, η)-PD). Let k ∈ N and η ∈ R≥0 . For
a training trace set S with |S| ≥ k, a synthetic trace y ∈R output by a generative model M with an input user
index d1 ∈ [|U|] is releasable with (k, η)-PD if there exist
at least k − 1 distinct training user indexes d2 , · · · , dk ∈[|U|]\{d1 } such that for any i, j ∈ [k],e−η p(y=M(dj )) ≤ p(y=M(di )) ≤ eη p(y=M(dj )). (1)The intuition behind (k, η)-PD can be explained as follows. Assume that user un is an input user of the synthetic trace y. Since y resembles the training trace snof un , it would be natural to consider an adversary who
attempts to recover sn (i.e., infer a pair of a user andthe whole sequence of events in sn ) from y. This attack
is called the tracking attack, and is decomposed into twophases: re-identification (or de-anonymization) and deobfuscation [62]. The adversary first uncovers the factthat user un is an input user of y, via re-identification.
Then she infers events of un via de-obfuscation. (k, η)PD can prevent re-identification because it guarantees
that the input user un is indistinguishable from at leastk − 1 other training users. Then the tracking attack isprevented even if de-obfuscation is perfectly done. A
large k and a small η are desirable for strong privacy.
(k, η)-PD can be used to alleviate the linkage of the
input user un and the synthetic trace y. However, y may
also leak information about parameters of the generative model Mn because y is generated using Mn . In
Section 3.5, we discuss the overall privacy of PPMTF
including this issue in detail.

5

Privacy-Preserving Multiple Tensor Factorization for Synthesizing Large-Scale Location Traces with Cluster-Specific Features

3.1 Overview
Proposed Method. Figure 2 shows an overview of
PPMTF (we formally define the symbols that newly appear in Figure 2 in Sections 3.2 to 3.4). It comprises thefollowing four steps.
(i). We compute a transition-count tensor RI and visitcount tensor RII from a training trace set S.
The transition-count tensor RI comprises the“User,” “Location,” and “Next Location” modes.
Its (n, i, j)-th element includes a transition-countof user un ∈ U from location xi ∈ X to xj ∈ X .
In other words, this tensor represents the movement pattern of each training user in the form oftransition-counts. The visit-count tensor RII comprises the “User,” “Location,” and “Time Slot”
modes. The (n, i, j)-th element includes a visitcount of user un at location xi in time slot lj ∈ L.
That is, this tensor includes a histogram of visited
locations for each user and each time slot.
(ii). We factorize the two tensors RI and RII simultaneously via MTF (Multiple Tensor Factorization)[35, 65], which factorizes multiple tensors into lowrank matrices called factor matrices along eachmode (axis). In MTF, one tensor shares a factor
matrix from the same mode with other tensors.
In our case, we factorize RI and RII into factor matrices A, B, C, and D, which respectively correspond to the “User,” “Location,” “Next Location,”and “Time Slot” mode. Here A and B are shared
between the two tensors. A, B, C, and D are parameters of our generative model, and therefore we callthem the MTF parameters. Let Θ = (A, B, C, D)
be the tuple of MTF parameters. We train MTF
parameters Θ from the two tensors via posterior
sampling [67], which samples Θ from its posterior
distribution given RI and RII .
(iii).We reconstruct two tensors from Θ. Then, given an
input user index n ∈ [|U|], we compute a transition-

er

(ii)

𝐑I

𝐑II

Next Location

Time Slot

(i)

𝐑I

𝐑II

Next Location Time Slot
𝐂
𝐃

(iii)
User unTime Transition-probability Visit-probability
Slot
Matrix
VectorQn,1𝝅n,11Training
Traces(v)
Synthetic
Traces𝐁

er

𝐀

Location Us

We propose PPMTF for synthesizing location traces. We
first present an overview (Section 3.1). Then we explain
the computation of two tensors (Section 3.2), the training of our generative model (Section 3.3), and the synthesis of traces (Section 3.4). Finally, we introduce the
PD (Plausible Deniability) test (Section 3.5).

MTF parameters 𝚯 = (𝐀, 𝐁, 𝐂, 𝐃)

Transition-count Visit-countTensorTensorLocation Us

3 Privacy-Preserving Multiple
Tensor Factorization (PPMTF)

(iv)

|ℒ|

Qn,|ℒ|

𝝅n,|ℒ|

Fig. 2. Overview of PPMTF with the following four steps: (i)computing a transition-count tensor and visit-count tensor, (ii)training MTF parameters via posterior sampling, (iii) computing
a transition-probability matrix and visit-probability vector via theMH algorithm and synthesizing traces, and (iv) the PD test.

probability matrix Qn,i and visit-probability vector πn,i of user un ∈ U for each time slot li ∈ L.
We compute them from the reconstructed tensors
via the MH (Metropolis-Hastings) algorithm [50],which modifies the transition matrix so that πn,i isa stationary distribution of Qn,i . Then we generate
a synthetic trace y ∈ R by using Qn,i and πn,i .
(iv).Finally, we perform the PD test [9], which verifies
whether y is releasable with (k, η)-PD.
We explain steps (i), (ii), (iii), and (iv) in Sections 3.2,
3.3, 3.4, and 3.5, respectively. We also explain how totune hyperparameters (parameters to control the training process) in PPMTF in Section 3.5. Below we explain
the utility, privacy, and scalability of PPMTF.
Utility. PPMTF achieves high utility by modeling statistical features of training traces using two tensors.
Specifically, the transition-count tensor represents the
movement pattern of each user in the form of transitioncounts, whereas the visit-count tensor includes a histogram of visited locations for each user and time slot.
Consequently, our synthetic traces preserve a timedependent population distribution, a transition matrix,
and a distribution of visit-counts per location; i.e., features (a), (b), and (c) in Section 1.
Furthermore, PPMTF automatically finds a clusterof users who have similar behaviors (e.g., those who always stay in Manhattan; those who often visit universities) and locations that are semantically similar (e.g.,
restaurants and bars) because factor matrices in tensor
factorization represent clusters [16]. Consequently, our
synthetic traces preserve the mobility behavior of sim-

Privacy-Preserving Multiple Tensor Factorization for Synthesizing Large-Scale Location Traces with Cluster-Specific Features

Privacy. PPMTF uses the PD test in [9] to provide PD
for a synthetic trace. In our experiments, we show that
PPMTF provides (k, η)-PD for reasonable k and η.
We also note that a posterior sampling-based
Bayesian learning algorithm, which produces a sample from a posterior distribution with bounded loglikelihood, provides DP without additional noise [67].
Based on this, we sample Θ from a posterior distribution given RI and RII to provide DP for Θ. However, theprivacy budget ε needs to be very large to achieve highutility in PPMTF. We discuss this issue in Appendix F.
Scalability. Finally, PPMTF achieves much higher
scalability than SGLT [8]. Specifically, the time complexity of [8] (semantic clustering) is O(|U|2 |X |3 |L|),which is very large for training traces with large |U|and |X |. On the other hand, the time complexity of
PPMTF is O(|U||X |2 ||L|) (see Appendix B for details),
which is much smaller than the synthesizer in [8]. In our
experiments, we evaluate the run time and show that
our method is applicable to much larger-scale training
datasets than SGLT.

3.2 Computation of Two Tensors
We next explain details of how to compute two tensors
from a training trace set S (i.e., step (i)).

se
r

u3

0 1 0
2 00 00 0
u1
x1 0 00 010 01 1
x2 1 00 001 20 0
x3 1 10 001 11 1
x4 1 10 00 0
x5 0 0 2
l1 l2 l3
Time SlotUU

0 0 1 0 0
0 10 00 00 00 0
u1
x1 0 00 000 000 000 00 2x2 0 00 100 010 000 00 0
x3 0 00 000 210 001 01 0
x4 0 00 00 00 10 0
x5 0 0 0 0 0x1 x2 x3 x4 x5
Next Locationu2

Visit-count tensor 𝐑II

u3

LocationserTransition-count tensor 𝐑I

Location

ilar users and the semantics of similar locations. They
also preserve a cluster-specific population distribution;i.e., feature (d) in Section 1,
More specifically, each column in A, B, C, andD represents a user cluster, location cluster, location
cluster, and time cluster, respectively. For example, elements with large values in the first column in B, C, andD may correspond to bars, bars, and night, respectively.
Then elements with large values in the first column in
A represent a cluster of users who go to bars at night.
In Section 4, we present visualization of some clusters, which can be divided into geographic clusters (e.g.,
north-eastern part of Tokyo) and semantic clusters (e.g.,trains, malls, universities). Semantic annotation of POIs
[19, 70] can also be used to automatically find what eachcluster represents (i.e., semantic annotation of clusters).
PPMTF also addresses sparseness of the tensors by
sharing A and B between the two tensors. It is shownin [65] that the utility is improved by sharing factor
matrices between tensors, especially when one of two
tensors is extremely sparse. We also confirmed that theutility is improved by sharing A and B.

6

u2

Fig. 3. Two tensors obtained from the training traces in Figure 1.

Two Tensors. Figure 3 presents an example of the two
tensors computed from the training traces in Figure 1.
The transition-count tensor includes a transition|U |×|X |×|X |count matrix for each user. Let RI ∈ Z≥0
be theItransition-count tensor, and rn,i,j∈ Z≥0 be its (n, i, j)I
th element. For example, r1,3,4= 2 in Figure 3 because
two transitions from x3 to x4 are observed in s1 of u1
in Figure 1. The visit-count tensor includes a histogramof visited locations for each user and each time slot. Let|U |×|X |×|L|II
RII ∈ Z≥0be the visit-count tensor, and rn,i,j∈IIZ≥0 be its (n, i, j)-th element. For example, r1,5,3 = 2
in Figure 3 because u1 visits x5 twice in l3 (i.e., fromtime instant 7 to 9) in Figure 1.
Let R = (RI , RII ). Typically, RI and RII are sparse;i.e., many elements are zeros. In particular, RI can be
extremely sparse because its size |RI | is quadratic in |X |.
Trimming. For both tensors, we randomly delete positive elements of users who have provided much morepositive elements than the average (i.e., outliers) in thesame way as [43]. This is called trimming, and is effective
for matrix completion [34]. The trimming is also used
to bound the log-likelihood in the posterior sampling
method [43] (we also show in Appendix F that the loglikelihood is bounded by the trimming). Similarly, we
set the maximum value of counts for each element, andtruncate counts that exceed the maximum number.
Specifically, let λI , λII ∈ N respectively represent the
maximum numbers of positive elements per user in RI
and RII . Typically, λI  |X | × |X | and λII  |X | × |L|.
For each user, if the number of positive elements in RI
exceeds λI , then we randomly select λI elements fromall positive elements, and delete the remaining positiveelements. Similarly, we randomly delete extra positiveIII∈ N be theelements in λII . In addition, let rmax
, rmaxmaximum counts for each element in RI and RII , reI
I
spectively. For each element, we truncate rn,i,j
to rmax
I
I
II
II
II
II
if rn,i,j > rmax (resp. rn,i,j to rmax if rn,i,j > rmax ).
In our experiments, we set λI = λII = 102 (as in [43])IIIand rmax= rmax= 10 because the number of positive

Privacy-Preserving Multiple Tensor Factorization for Synthesizing Large-Scale Location Traces with Cluster-Specific Featureselements per user and the value of counts were respectively less than 100 and 10 in most cases. In other words,
the utility does not change much by increasing the valIII . We also confirmed thatues of λI , λII , rmax
, and rmax
IIImuch smaller values (e.g., λI = λII = rmax
= rmax= 1)
result in a significant loss of utility.

7

𝜇0 , 𝛽0 , 𝑊0 , 𝜈0
𝚿𝐀

𝚿𝐁

𝚿𝐂

𝚿𝐃

𝐀

𝐁

𝐂

𝐃

𝐑𝐈

𝐑𝐈𝐈

MTF parameters 𝚯
Tensors 𝐑

𝛼

3.3 Training MTF Parameters
After computing R = (RI , RII ), we train the MTF parameters Θ = (A, B, C, D) via posterior sampling (i.e.,
step (ii)). Below we describe our MTF model and thetraining of Θ.
Model. Let z ∈ N be the number of columns (factors)
in each factor matrix. Let A ∈ R|U |×z , B ∈ R|X |×z , C ∈R|X |×z , and D ∈ R|L|×z be the factor matrices. Typically, the number of columns is much smaller than thenumbers of users and locations; i.e., z  min{|U|, |X |}.
In our experiments, we set z = 16 as in [49] (we alsochanged the number z of factors from 16 to 32 and confirmed that the utility was not changed much).
Let ai,k , bi,k , ci,k , di,k ∈ R be the (i, k)-th elementsof A, B, C, and D, respectively. In addition, let R̂I ∈
R|U |×|X |×|X | and R̂II ∈ R|U |×|X |×|L| respectively represent two tensors that can be reconstructed from Θ.
IIISpecifically, let r̂n,i,j∈ R and r̂n,i,j∈ R be the (n, i, j)III
th elements of R̂ and R̂ , respectively. Then R̂I andR̂II are given by:X
X
I
IIr̂n,i,j=an,k bi,k cj,k , r̂n,i,j=an,k bi,k dj,k , (2)
k∈[z]

k∈[z]

where A and B are shared between R̂I and R̂II .
For MTF parameters Θ, we use a hierarchical Bayesmodel [58] because it outperforms the non-hierarchicalone [57] in terms of the model’s predictive accuracy.
Specifically, we use a hierarchical Bayes model shown in
Figure 4. Below we explain this model in detail.
For the conditional distribution p(R|Θ) of the two
tensors R = (RI , RII ) given the MTF parametersIΘ = (A, B, C, D), we assume that each element rn,i,j
II ) is independently generated from a normal(resp. rn,i,jIII ) and precision
distribution with mean r̂n,i,j(resp. r̂n,i,j(reciprocal of the variance) α ∈ R≥0 . In our experiments, we set α to various values from 10−6 to 103 .
Here we randomly select a small number of zero elements in R to improve the scalability in the same way
as [3, 54]. Specifically, we randomly select ρI ∈ N andρII ∈ N zero elements for each user in RI and RII , respectively, where ρI  |X | × |X | and ρII  |X | × |L|(in our experiments, we set ρI = ρII = 103 ). We

Fig. 4. Graphical model of PPMTF.

I
treat the remaining zero elements as missing. Let In,i,j
II ) be the indicator function that takes 0 if
(resp. In,i,jIII ) is missing, and takes 1 otherwise.
rn,i,j (resp. rn,i,jIII ) takes 1 at most λI + ρI
Note that In,i,j (resp. In,i,jIIII(resp. λ + ρ ) elements for each user, where λI (resp.
λII ) is the maximum number of positive elements peruser in RI (resp. RII ).
Then the distribution p(R|Θ) can be written as:p(R|Θ) = p(RI |A, B, C)p(RII |A, B, D)
Y
I
I
I
=
[N (rn,i,j|r̂n,i,j, α−1 )]In,i,j
n,i,j

·

Y

II

II
II
[N (rn,i,j|r̂n,i,j, α−1 )]In,i,j ,

(3)

n,i,j

where N (r|µ, α−1 ) denotes the probability of r in thenormal distribution with mean µ and precision α (i.e.,variance α−1 ).
Let ai , bi , ci , di ∈ Rz be the i-th rows of A, B, C,
and D, respectively. For a distribution of Θ = (A, B,
C, D), we assume the multivariate normal distribution:Q
p(A|ΨA ) = n N (an |µA , Λ−1A)Q−1p(B|ΨB ) = n N (bn |µB , ΛB )Q
p(C|ΨC ) = n N (cn |µC , Λ−1C )Qp(D|ΨD ) = n N (dn |µD , Λ−1D ),where µA , µB , µC , µD ∈ Rz are mean vectors, ΛA ,
ΛB , ΛC , ΛD ∈ Rz×z are precision matrices, and ΨA =(µA , ΛA ), ΨB = (µB , ΛB ), ΨC = (µC , ΛC ), ΨD =(µD , ΛD ).
The hierarchical Bayes model assumes a distribution for each of ΨA , ΨB , ΨC , and ΨD , which is called a
hyperprior. We assume ΨZ ∈ {ΨA , ΨB , ΨC , ΨD } follows
a normal-Wishart distribution [10], i.e., the conjugate
prior of a multivariate normal distribution:p(ΨZ ) = p(µZ |ΛZ )p(ΛZ )= N (µZ |µ0 , (β0 ΛZ )−1 )W(ΛZ |W0 , ν0 ),(4)

where µ0 ∈ Rz , β0 ∈ R, and W(Λ|W0 , ν0 ) denotes theprobability of Λ ∈ Rz×z in the Wishart distribution

8

Privacy-Preserving Multiple Tensor Factorization for Synthesizing Large-Scale Location Traces with Cluster-Specific Features

Posterior Sampling of Θ. We train Θ based on theposterior sampling method [67]. This method trains Θfrom R by sampling Θ from the posterior distributionp(Θ|R). To sample Θ from p(Θ|R), we use Gibbs sampling [50], which samples each variable in turn, conditioned on the current values of the other variables.
Specifically, we sample ΨA , ΨB , ΨC , ΨD , A, B, C,
and D in turn. We add superscript “(t)” to these variables to denote the sampled values at the t-th iteration.
For initial values with “(0)”, we use a random initialization method [4] that initializes each element as a random number in [0, 1] because it is widely used. Then,
(t)
(t)
(t)
(t)
we sample ΨA , ΨB , ΨC , ΨD , A(t) , B(t) , C(t) , andD(t) from the conditional distribution given the current
values of the other variables, and iterate the sampling
for a fixed number of times (we omit the details of thesampling algorithm for lack of space).
Gibbs sampling guarantees that the sampling distributions of A(t) , · · · , D(t) approach the posterior distributions p(A|R), · · · , p(D|R) as t increases. Therefore,
Θ(t) = (A(t) , B(t) , C(t) , D(t) ) approximates Θ sampledfrom the posterior distribution p(Θ|R) for large t. In
our experiments, we discarded the first 99 samples as“burn-in”, and used Θ(100) as an approximation of Θ.
We also confirmed that the model’s predictive accuracy
converged within 100 iterations.

3.4 Generating Traces via MH
After training Θ = (A, B, C, D), we generate synthetic
traces via the MH (Metropolis-Hastings) algorithm [50](i.e., step (iii)). Specifically, given an input user indexn ∈ [|U|], we generate a synthetic trace y ∈ R that
resembles sn of user un ∈ U from (an , B, C, D). In otherwords, the parameters of the generative model Mn ofuser un are (an , B, C, D).
Let Q be the set of |X | × |X | transition-probability
matrices, and C be the set of |X |-dimensional probability
vectors (i.e., probability simplex). Given a transitionprobability matrix Q ∈ Q and a probability vector π ∈C, the MH algorithm modifies Q to Q0 ∈ Q so thatthe stationary distribution of Q0 is equal to π. Q is aconditional distribution called a proposal distribution,and π is called a target distribution.

𝐐𝑛 ,1

𝐑I𝒏

Next Location

𝝅n,1

𝐐∗𝑛
𝐐𝑛 ,|ℒ|

𝝅n,|ℒ|

LocationTransition-prob. Visit-prob.
Matrix
VectorLocationwith parameters W0 ∈ Rz×z and ν0 ∈ R (W0 and ν0
represent the scale matrix and the number of degrees of
freedom, respectively). µ0 , β0 , W0 , and ν0 are parameters of the hyperpriors, and are determined in advance.
In our experiments, we set µ0 = 0, β0 = 2, ν0 = z, andW0 to the identity matrix, in the same way as [58].

𝐑II𝒏Time Slot

Fig. 5. Computation of (Qn,i , πn,i ) via MH. We compute Q∗n
I , and πIIfrom R̂nn,i from R̂n . Then for each time slot li ∈ L, we
modify Q∗n to Qn,i whose stationary distribution is πn,i .

In step (iii), given the input user index n ∈ [|U|], wereconstruct the transition-count matrix and visit-countmatrix of user un , and use the MH algorithm to makea transition-probability matrix of un consistent with avisit-probability vector of un for each time slot. Figure 5
shows its overview. Specifically, let R̂nI ∈ R|X |×|X | andR̂nII ∈ R|X |×|L| be the n-th matrices in R̂I and R̂II ,respectively (i.e., reconstructed transition-count matrix
and visit-count matrix of user un ). We first compute R̂nI
and R̂nII from (an , B, C, D) by (2). Then we compute atransition-probability matrix Q∗n ∈ Q of user un fromR̂nI by normalizing counts to probabilities. Similarly, we
compute a visit-probability vector πn,i ∈ C of user unfor each time slot li ∈ L from R̂nII by normalizing counts
to probabilities. Then, for each time slot li ∈ L, we
modify Q∗n to Qn,i ∈ Q via the MH algorithm so thatthe stationary distribution of Qn,i is equal to πn,i . Thenwe generate a synthetic trace using (Qn,i , πn,i ).
Below we explain step (iii) in more detail.
Computing (Qn,i , πn,i ) via MH. We first computethe n-th matrix R̂nI ∈ R|X |×|X | in R̂I from Θ by (2).
Then we compute Q∗n ∈ Q from R̂nI by normalizing
counts to probabilities as explained below. We assign
a very small positive value φ ∈ R≥0 (φ = 10−8 inour experiments) to elements in R̂nI with values smaller
than φ. Then we normalize R̂nI to Q∗n so that the sum
over each row in Q∗n is 1. Since we assign φ (= 10−8 )to elements with smaller values in R̂nI , the transitionprobability matrix Q∗n is regular [50]; i.e., it is possible
to get from any location to any location in one step. This
allows πn,i to be the stationary distribution of Qn,i , as
explained later in detail.
We then compute the n-th matrix R̂nII ∈ R|X |×|L| in
R̂II from Θ by (2). For each time slot li ∈ L, we assign φ(= 10−8 ) to elements with smaller values in R̂nII . Then
we normalize the i-th column of R̂nII to πn,i ∈ C so thatthe sum of πn,i is one.
We use Q∗n as a proposal distribution and πn,i asa target distribution, and apply the MH algorithm to9Privacy-Preserving Multiple Tensor Factorization for Synthesizing Large-Scale Location Traces with Cluster-Specific Features

User
un

y

x3

x2

Synthetic Tracex4 x1 x1 x2

Distribution 𝝅n,1𝐐𝑛 ,1Hour (AM) 7
Time Instant 1
Time Slot

2

3

x3

𝐐𝑛 ,2
8
4

l1

5

l2

x5

x5

𝐐𝑛 ,3
6

9
7

8

9

l3

Fig. 6. Generation of a synthetic trace (|X | = 5, |T | = 9, |L| =3). Each location is randomly generated from a distribution in thesame time slot.

obtain a transition-probability matrix Qn,i whose stationary distribution is πn,i . For Q ∈ Q and a, b ∈ [|X |],we denote by Q(xb |xa ) ∈ [0, 1] the transition probability
from xa ∈ X to xb ∈ X (i.e., the (a, b)-th element of Q).
Similarly, given π ∈ C, we denote by π(xa ) ∈ [0, 1] thevisit probability at xa ∈ X . Then, the MH algorithm
computes Qn,i (xb |xa ) for xa 6= xb as follows:
πn,i (xb )Q∗n (xa |xb )Qn,i (xb |xa ) = Q∗n (xb |xa ) min 1, πn,i(xa )Q∗ (xb |xa ) ,n(5)
and computes Qn,i (xa |xa ) as follows: Qn,i (xa |xa ) =P
1 − b6=a Qn,i (xb |xa ). Note that Qn,i is regular because
all elements in Q∗n and πn,i are positive. Then the MH
algorithm guarantees that πn,i is a stationary distribution of Qn,i [50].
Generating Traces. After computing (Qn,i , πn,i ) via
the MH algorithm, we synthesize a trace y ∈ R of user
un as follows. We randomly generate the first location intime slot l1 from the visit-probability distribution πn,1 .
Then we randomly generate the subsequent location intime slot li using the transition-probability matrix Qn,i .
Figure 6 shows an example of synthesizing a trace y ofuser un . In this example, a location at time instant 7
is randomly generated from the conditional distribution
Qn,3 given the location x2 at time instant 6.
The synthetic trace y is generated in such a waythat a visit probability in time slot li is given by πn,i .
In addition, the transition matrix is computed by using Q∗n as a proposal distribution. Therefore, we cansynthesize traces that preserve the statistical feature oftraining traces such as the time-dependent population
distribution and the transition matrix.

3.5 Privacy Protection
We finally perform the PD test for a synthetic trace y.
Let MPPMTF be our generative model in step (iii)that, given an input user index n ∈ [|U|], outputs a syn-thetic trace y ∈ R with probability p(y = MPPMTF (n)).
Let σ : T → X be a function that, given time instantt ∈ T , outputs an index of the location at time instantt in y; e.g., σ(1) = 3, σ(2) = 2, · · · , σ(9) = 5 in Figure 6.
Furthermore, let ω : T → L be a function that, given
time instant t ∈ T , outputs an index of the corresponding time slot; e.g., ω(1) = ω(2) = ω(3) = 1, · · · , ω(7) =ω(8) = ω(9) = 3 in Figure 6.
Recall that the first location in y is randomly generated from πn,1 , and the subsequent location at time
instant t ∈ T is randomly generated from Qn,ω(t) . Then,p(y = MPPMTF (n))
Q|T |= πn,1 (xσ(1) ) t=2 Qn,ω(t) (xσ(t) |xσ(t−1) ).
Thus, given y ∈ R, we can compute p(y =MPPMTF (m)) for any m ∈ [|U|] as follows: (i) compute (Qm,i , πm,i ) for each time slot li ∈ L via the MH
algorithm (as described in Section 3.4); (ii) computep(y = MPPMTF (m)) using (Qm,i , πm,i ). Then we canverify whether y is releasable with (k, η)-PD by counting the number of training users such that (1) holds.
Specifically, we use the following PD test in [9]:Privacy Test 1 (Deterministic Test in [9]). Let k ∈ Nand η ∈ R≥0 . Given a generative model M, training
user set U, input user index n ∈ [|U|], and synthetic
trace y, output pass or fail as follows:
1. Let i ∈ Z≥0 be a non-negative integer that satisfies:e−(i+1)η < p(y = M(n)) ≤ e−i η .

(6)

2. Let k 0 ∈ Z≥0 be the number of training user indexesm ∈ [|U|] such that:e−(i+1)η < p(y = M(m)) ≤ e−i η .

(7)

3. If k 0 ≥ k, then return pass, otherwise return fail.
By (1), (6), and (7), if y passes Privacy Test 1, theny is releasable with (k, η)-PD. In addition, (k, η)-PD isguaranteed even if Θ is not sampled from the exact posterior distribution p(Θ|R).
The time complexity of Privacy Test 1 is linear in|U|. In this paper, we randomly select a subset U ∗ ⊆ U
of training users from U (as in [9]) to ascertain more
quickly whether k 0 ≥ k or not. Specifically, we initialize
k 0 to 0, and check (7) for each training user in U ∗ ∪ {un }(increment k 0 if (7) holds). If k 0 ≥ k, then we return
pass (otherwise, return fail). The time complexity ofthis faster version of Privacy Test 1 is linear in |U ∗ |(≤ |U|). A smaller |U ∗ | leads to a faster (k, η)-PD test at
the expense of fewer synthetic traces passing the test.

Privacy-Preserving Multiple Tensor Factorization for Synthesizing Large-Scale Location Traces with Cluster-Specific Features

In Section 4, we use the faster version of Privacy
Test 1 with |U ∗ | = 32000, k = 10 to 200, and η = 1 toguarantee (k, η)-PD for reasonable k and η (note thatε = 1 is considered to be reasonable in ε-DP [23, 38]).
Overall Privacy. As described in Section 2.2, even if
a synthetic trace y satisfies (k, η)-PD, y may leak information about the MTF parameters. We finally discuss
the overall privacy of y including this issue.
Given the input user index n, PPMTF generates yfrom (an , B, C, D), as described in Section 3.4. Since thelinkage of the input user un and y is alleviated by PD,
the leakage of an is also alleviated by PD. Therefore,
the remaining issue is the leakage of (B, C, D).
Here we note that B and C are information aboutlocations (i.e., location profiles), and D is information
about time (i.e., time profile). Thus, even if the adversary perfectly infers (B, C, D) from y, it is hard to inferprivate information (i.e., training traces S) of users U
from (B, C, D) (unless she obtains user profile A). In
fact, some studies on privacy-preserving matrix factorization [45, 53] release an item profile publicly. Similarly, SGLT [8] assumes that semantic clusters of locations (parameters of their generative model) leak almost
no information about U because the location clusters area kind of location profile. We also assume that the location and time profiles leak almost no information about
users U. Further analysis is left for future work.
Tuning Hyperparameters. As described in SecI
tions 3.2, 3.3, and 3.5, we set λI = λII = 102 , rmax=
II
rmax = 10 (because the number of positive elements
per user and the value of counts were respectively lessthan 100 and 10 in most cases), z = 16 (as in [49]),
ρI = ρII = 103 , |U ∗ | = 32000, and changed α from 10−6
and 103 in our experiments. If we set these values to verysmall values, the utility is lost (we show its example bychanging α in our experiments). For the parameters of
the hyperpriors, we set µ0 = 0, β0 = 2, ν0 = z, and W0
to the identity matrix in the same way as [58].
We set the hyperparameters as above based on theprevious work or the datasets. To optimize the hyperparameters, we could use, for example, cross-validation[10], which assesses the hyperparameters by dividing adataset into a training set and testing (validation) set.

4 Experimental Evaluation
In our experiments, we used two publicly available
datasets: the SNS-based people flow data [52] and theFoursquare dataset in [68]. The former is a relatively

10small-scale dataset with no missing events. It is used tocompare the proposed method with two state-of-the-art
synthesizers [8, 9]. The latter is one of the largest publicly available location datasets; e.g., much larger than
[14, 55, 69, 74]. Since the location synthesizer in [8] cannot be applied to this large-scale dataset (as shown in
Section 4.4), we compare the proposed method with [9].

4.1 Datasets
SNS-Based People Flow Data. The SNS-based people flow data [52] (denoted by PF) includes artificial
traces around the Tokyo metropolitan area. The traces
were generated from real geo-tagged tweets by interpolating locations every five minutes using railway androad information [59].
We divided the Tokyo metropolitan area into 20×20
regions; i.e., |X | = 400. Then we set the interval between
two time instants to 20 minutes, and extracted tracesfrom 9:00 to 19:00 for 1000 users (each user has a single
trace comprising 30 events). We also set time slots to
20 minutes long from 9:00 to 19:00. In other words, we
assumed that each time slot comprises one time instant;i.e., |L| = 30. We randomly divided the 1000 traces into
500 training traces and 500 testing traces; i.e., |U| = 500.
The training traces were used for training generativemodels and synthesizing traces. The testing traces wereused for evaluating the utility.
Since the number of users is small in PF, we generated ten synthetic traces from each training trace (eachsynthetic trace is from 9:00 to 19:00) and averaged theutility and privacy results over the ten traces to stabilize
the performance.
Foursquare Dataset.
The Foursquare dataset(Global-scale Check-in Dataset with User Social Networks) [68] (denoted by FS) includes 90048627 real
check-ins by 2733324 users all over the world.
We selected six cities with numerous check-ins andwith cultural diversity in the same way as [68]: Istanbul (IST), Jakarta (JK), New York City (NYC), Kuala
Lumpur (KL), San Paulo (SP), and Tokyo (TKY). For
each city, we extracted 1000 POIs, for which the numberof visits from all users was the largest; i.e., |X | = 1000.
We set the interval between two time instants to 1 hour
(we rounded down minutes), and assigned every 2 hours
into one of 12 time slots l1 (0-2h), · · · , l12 (22-24h) ina cyclic manner; i.e., |L| = 12. For each city, we randomly selected 80% of traces as training traces and used
the remaining traces as testing traces. The numbers|U| of users in IST, JK, NYC, KL, SP, and TKY werePrivacy-Preserving Multiple Tensor Factorization for Synthesizing Large-Scale Location Traces with Cluster-Specific Features

219793, 83325, 52432, 51189, 42100, and 32056, respectively. Note that there were many missing events in FS
because FS is a location check-in dataset. The numbersof temporally-continuous events in the training tracesof IST, JK, NYC, KL, SP, and TKY were 109027, 19592,
7471, 25563, 13151, and 47956, respectively.
From each training trace, we generated one synthetic trace with the length of one day.

4.2 Location Synthesizers
We evaluated the proposed method (PPMTF), the synthetic location traces generator in [8] (SGLT), and thesynthetic data generator in [9] (SGD).
III=In PPMTF, we set λI = λII = 102 , rmax
= rmaxIII
310, z = 16, ρ = ρ = 10 , µ0 = 0, β0 = 2, ν0 = z, andW0 to the identity matrix, as explained in Section 3.
Then we evaluated the utility and privacy for each value.
In SGLT [8], we used the SGLT tool (C++) in [7].
We set the location-removal probability parc to 0.25, thelocation merging probability parm to 0.75, and the randomization multiplication factor parv to 4 in the same
way as [8] (for details of the parameters in SGLT, see [8]).
For the number c of semantic clusters, we attempted
various values: c = 50, 100, 150, or 200 (as shown later,
SGLT provided the best performance when c = 50 or100). For each case, we set the probability parl of removing the true location in the input user to various
values from 0 to 1 (parl = 1 in [8]) to evaluate thetrade-off between utility and privacy.
In SGD [9], we trained the transition matrix for eachtime slot (|L| × |X | × |X | elements in total) and thevisit-probability vector for the first time instant (|X | elements in total) from the training traces via maximum
likelihood estimation. Note that the transition matrix
and the visit-probability vector are common to all users.
Then we generated a synthetic trace from an input user
by copying the first ξ ∈ Z≥0 events of the input user andgenerating the remaining events using the trained transition matrix. When ξ = 0, we randomly generated a location at the first time instant using the visit-probability
vector. For more details of SGD for location traces, see
Appendix C. We implemented PPMTF and SGD with
C++, and published it as open-source software [1].

4.3 Performance Metrics
Utility. In our experiments, we evaluated the utility
listed in Section 1.

11(a) Time-Dependent Population Distribution. We
computed a frequency distribution (|X |-dim vector) of
the testing traces and that of the synthetic traces foreach time slot. Then we evaluated the average total variation between the two distributions over all time slots(denoted by TP-TV).
Frequently visited locations are especially important for some tasks [19, 75]. Therefore, for each timeslot, we also selected the top 50 locations, whose frequencies in the testing traces were the largest, and regarded the absolute error for the remaining locations in
TP-TV as 0 (TP-TV-Top50).
(b) Transition Matrix. We computed an averagetransition-probability matrix (|X | × |X | matrix) over allusers and all time instances from the testing traces. Similarly, we computed an average transition-probability
matrix from the synthetic traces.
Since each row of the transition matrix represents a
conditional distribution, we evaluated the EMD (EarthMover’s Distance) between the two conditional distributions over the x-axis (longitude) and y-axis (latitude),and averaged it over all rows (TM-EMD-X and TMEMD-Y). TM-EMD-X and TM-EMD-Y represent how
the two transition matrices differ over the x-axis andy-axis, respectively. They are large especially when one
matrix allows only a transition between close locations
and the other allows a transition between far-away locations (e.g., two countries). The EMD is also used in[8] to measure the difference in two transition matrices.
We did not evaluate the two-dimensional EMD, because
the computational cost of the EMD is expensive.
(c) Distribution of Visit-Fractions. Since we used
POIs in FS (regions in PF), we evaluated how well thesynthetic traces preserve a distribution of visit-fractions
in FS. We first excluded testing traces that have a few
events (fewer than 5). Then, for each of the remaining traces, we computed a fraction of visits for each
POI. Based on this, we computed a distribution of visitfractions for each POI by dividing the fraction into 241
21], ( 24, 24], · · · , ( 23bins as (0, 24
24 , 1). Similarly, we computed a distribution of visit-fractions for each POI fromthe synthetic traces. Finally, we evaluated the total variation between the two distributions (VF-TV).
(d) Cluster-Specific Population Distribution. To
show that PPMTF is also effective in this respect, we
conducted the following analysis. We used the fact thateach column in the factor matrix A represents a cluster(z = 16 clusters in total). Specifically, for each column
in A, we extracted the top 10% users whose values in
the column are the largest. These users form a clus-

Privacy-Preserving Multiple Tensor Factorization for Synthesizing Large-Scale Location Traces with Cluster-Specific Features

ter who exhibit similar behavior. For some clusters, we
visualized factor matrices and the frequency distributions (i.e., cluster-specific population distributions) of
the training traces and synthetic traces.
Privacy. In PF, we evaluated the three synthesizers.
Although PPMTF and SGD provide (k, η)-PD in Definition 1, SGLT provides PD using a semantic distance
between traces [8], which differs from PD in Definition 1.
To compare the three synthesizers using the same
privacy metrics, we considered two privacy attacks: reidentification (or de-anonymization) attack [26, 47, 62]
and membership inference attack [31, 61]. In the reidentification attack, the adversary identifies, for eachsynthetic trace y, an input user of y from |U| = 500
training users. We evaluated a re-identification rate as
the proportion of correctly identified synthetic traces.
In the membership inference attack, the adversary
obtains all synthetic traces. Then the adversary determines, for each of 1000 users (500 training users and
500 testing users), whether her trace is used for training
the model. Here training users are members and testingusers are non-members (they are randomly chosen, as
described in Section 4.1). We used membership advantage [71] as a privacy metric in the same way as [31].
Specifically, let tp, tn, f p, and f n be the number oftrue positives, true negatives, false positives, and false
negatives, respectively, where “positive/negative” represents a member/non-member. Then membership advantage is defined in [71] as the difference between the true
positive rate and the false positive rate; i.e., membertp
fp
tp−f pship advantage = tp+f
n − f p+tn =500 . Note thatmembership advantage can be easily translated into
membership inference accuracy, which is the proportp+tn
tion of correct adversary’s outputs (= tp+tn+f
p+f n =
tp+tn
1000 ), as follows: membership inference accuracy =membership advantage+1(since tn + f p = 500). A random
2

guess that randomly outputs “member” with probability q ∈ [0, 1] achieves advantage = 0 and accuracy = 0.5.
For both the re-identification attack and membership inference attack, we assume the worst-case scenario
about the background knowledge of the adversary; i.e.,
maximum-knowledge attacker model [20]. Specifically,
we assumed that the adversary obtains the 1000 original
traces (500 training traces and 500 testing traces) in PF.
Note that the adversary does not know which ones aretraining traces (and therefore performs the membership
inference attack). The adversary uses the 1000 original
traces to build an attack model. For a re-identification
algorithm, we used the Bayesian re-identification algorithm in [47]. For a membership inference algorithm, we

12

implemented a likelihood-ratio based membership inference algorithm, which partly uses the algorithm in [48].
For details of the attack algorithms, see Appendix D.
Note that evaluation might be difficult for a partialknowledge attacker who has less background knowledge.
In particular, when the amount of training data is small,
it is very challenging to accurately train an attack model(transition matrices) [46–48]. We note, however, that if
a location synthesizer is secure against the maximumknowledge attacker, then we can say that it is also secureagainst the partial-knowledge attacker, without implementing clever attack algorithms. Therefore, we focus
on the maximum-knowledge attacker model.
In FS, we used (k, η)-PD in Definition 1 as a privacy
metric because we evaluated only PPMTF and SGD. As
a PD test, we used the (faster) Privacy Test 1 with|U ∗ | = 32000, k = 10 to 200, and η = 1.
Scalability. We measured the time to synthesize traces
using the ABCI (AI Bridging Cloud Infrastructure) [51],which is a supercomputer ranking 8th in the Top 500(as of June 2019). We used one computing node, which
consists of two Intel Xeon Gold 6148 processors (2.40
GHz, 20 Cores) and 412 GB main memory.

4.4 Experimental Results in PF
Utility and Privacy.
Figure 7 shows the reidentification rate, membership advantage, and utility
with regard to (a) the time-dependent population distribution and (b) transition matrix in PF. Here, we set the
precision α in PPMTF to various values from 0.5 to 1000.
Uniform represents the utility when all locations in synthetic traces are independently sampled from a uniform
distribution. Training represents the utility of the training traces; i.e., the utility when we output the training
traces as synthetic traces without modification. Ideally,
the utility of the synthetic traces should be much better
than that of Uniform and close to that of Training.
Figure 7 shows that PPMTF achieves TP-TV and
TP-TV-Top50 close to Training for while protecting user
privacy. For example, PPMTF achieves TP-TV = 0.43
and TP-TV-Top50 = 0.13, both of which are close to
those of Training (TP-TV = 0.39 and TP-TV-Top50 =0.12), while keeping re-identification rate < 0.02 andmembership advantage < 0.055 (membership inference
accuracy < 0.53). We consider that PPMTF achieved
low membership advantage because (1) held for not onlyk − 1 training users but testing users (non-members).
In SGLT and SGD, privacy rapidly gets worse withdecrease in TP-TV and TP-TV-Top50. This is because

13

Privacy-Preserving Multiple Tensor Factorization for Synthesizing Large-Scale Location Traces with Cluster-Specific Features

0.4

0.6

0.8

1

0.125
0.1

0

0.2

0.4

0.6

0.8

Re-identification Rate

1

0

7
6
5
4
3
2
1
00

0.2

0.4

0.6

0.8

Re-identification Rate

1

0.45
0.4
0.35

0.2

0.4

0.6

0.8

Re-identification Rate

1

7
6
5
4
3
2
1
0

Uniform
Training

0.2

TP-TV-Top50

0.15

SGLT (150)
SGLT (200)

SGLT (50)
SGLT (100)
0.175

0

0.2

0.4

0.6

0.8

Membership Advantage

1

0.15
0.125
0.1

TM-EMD-Y

0.2

Re-identification Rate

TP-TV

0.175

PPMTF
SGD

0.65
0.6
0.55
0.5

TM-EMD-X

7
6
5
4
3
2
1
0

Uniform
Training

0.2

TM-EMD-Y

0.45
0.4
0.35
0

(ii) Membership Inference

SGLT (150)
SGLT (200)

SGLT (50)
SGLT (100)TP-TV-Top50

0.65
0.6
0.55
0.5

TM-EMD-X

TP-TV(i) Re-identification
PPMTF
SGD

0

0.2

0.4

0.6

0.8

Membership Advantage

1

7
6
5
4
3
2
1
0

0

0

0.2

0.4

0.6

0.8

1

0.2

0.4

0.6

0.8

1

Membership Advantage

Membership Advantage

Fig. 7. Privacy and utility in PF. The number in SGLT represents the number c of clusters. In PPMTF, SGLT and SGD, we varied α,parl and ξ, respectively. Lower is better in all of the utility metrics.

both SGLT and SGD synthesize traces by copying oversome events from the training traces. Specifically, SGLT
(resp. SGD) increases the number of copied events bydecreasing parl (resp. increasing ξ). Although a larger
number of copied events result in a decrease of both
TP-TV and TP-TV-Top50, they also result in the rapid
increase of the re-identification rate. This result is consistent with the uniqueness of location data; e.g., only
three locations are sufficient to uniquely characterize
80% of the individuals among 1.5 million people [18].
Figure 7 also shows that PPMTF performs worse
than SGLT and SGD in terms of TM-EMD-X and TMEMD-Y. This is because PPMTF modifies the transition
matrix so that it is consistent with a visit-probability
vector using the MH algorithm (SGLT and SGD do
not modify the transition matrix). It should be noted,
however, that PPMTF significantly outperforms Uniformwith regard to TM-EMD-X and TM-EMD-Y. This meansthat PPMTF preserves the transition matrix well.
Analysis on Cluster-Specific Features. Next, weshow the utility with regard to (d) the cluster-specificpopulation distribution. Specifically, we show in Figure 8 the frequency distributions of training traces andsynthetic traces and the columns of factor matrices B
and C for three clusters (we set α = 200 because it provided almost the best utility in Figure 7; we also normalized elements in each column of B and C so thatthe square-sum is one). Recall that for each cluster, we
extracted the top 10% users; i.e., 50 users.
Figure 8 shows that the frequency distributions of
training traces differ from cluster to cluster, and that

Training traces Synthetic traces Factor matrix B Factor matrix C

0.010
0.008

(i)

0.006
0.004
0.002
0.000
0.010
0.008

(ii)

0.006
0.004
0.002
0.000
0.010
0.008

(iii)

0.006
0.004
0.002
0.000

Fig. 8. Frequency distributions and the columns of factor matrices B and C for three clusters (50 users for each cluster) in
PF. The green line in (ii) and the orange line in (iii) represent
subways (Shinjuku and Fukutoshin lines, respectively).

the users in each cluster exhibit similar behavior; e.g.,
the users in (i) stay in the northeastern area of Tokyo;
the users in (ii) and (iii) often use the subways. PPMTF
models such a cluster-specific behavior via B and C,
and synthesizes traces that preserve the behavior using
B and C. Figure 8 shows that PPMTF is useful for geodata analysis such as modeling human location patterns[40] and map inference [5, 41].
Scalability. We also measured the time to synthesize traces from training traces. Here we generated one
synthetic trace from each training trace (500 synthetic
traces in total), and measured the time. We also changed
the numbers of users and locations (i.e., |U|, |X |) for various values from 100 to 1000 to see how the running timedepends on |U| and |X |.

Privacy-Preserving Multiple Tensor Factorization for Synthesizing Large-Scale Location Traces with Cluster-Specific Features

TP-TV

4.5 Experimental Results in FS
Utility and Privacy. In FS, we set α = 200 in PPMTF(as in Figures 8 and 9). In SGD, we set ξ = 0 for thefollowing two reasons: (1) the re-identification rate ishigh for ξ ≥ 1 in Figure 7 because of the uniqueness oflocation data [18]; (2) the event in the first time slot

w/o PDk = 200

k = 10Training

IST

NYC

TKY
TP-TV-Top50

0.2
0.15
0.1
0.05
0

IST

NYC

IST

IST JK NYC KLSP TKYIST JK NYC KLSP TKYIST JK NYC KLSP TKYIST JK NYC KLSP TKYIST JK NYC KLSP TKY

0.1
0.05

TKY

0
5

NYC

TKY

4
3
2
1
0
5

IST

NYC

TKY

4
3
2
1
0
0.75

0.75

VF-TV

0.6
0.45
0.3

0.6
0.45
0.3
0.15

0.15
0

SGD
Training

0.2

TM-EMD-X

3.5
3
2.5
2
1.5
1
0.5
0

Uniform
PPMTF (k = 10)

0.15

TM-EMD-Y

2.5
2
1.5
1
0.5
0

0.6
0.5
0.4
0.3
0.2
0.1
0

0.25

3

TM-EMD-X

Figure 9 shows the results (we set α = 200 in
PPMTF, and c = 100 and parl = 1 in SGLT; we also obtained almost the same results for other values). Here we
excluded the running time of SGD because it was verysmall; e.g., less than one second when |U| = 1000 and|X | = 400 (we compare the running time of PPMTF withthat of SGD in FS, as described later). The running timeof SGLT is much larger than that of PPMTF. Specifically, the running time of SGLT is quadratic in |U| (e.g.,when |X | = 400, SGLT(T) requires 0.47 and 47 hours
for |U| = 100 and 1000, respectively) and cubic in |X |(e.g., when |U| = 100, SGLT(T) requires 8.1 × 10−3 and8.4 hours for |X | = 100 and 1000, respectively). On theother hand, the running time of PPMTF is linear in |U|(e.g., PPMTF(S) requires 6.3×10−5 and 5.9×10−4 hours
for |U| = 100 and 1000, respectively) and quadratic in|X | (e.g., PPMTF(S) requires 9.3 × 10−3 and 0.96 hours
for |X | = 100 and 1000, respectively). This is consistent
with the time complexity described in Section 3.1.
From Figure 9, we can estimate the running timeof SGLT for generating large-scale traces. Specifically,
when |U| = 219793 and |X | = 1000 as in IST of
FS, SGLT(T) (semantic clustering) would require about4632 years (=8.4 × (219793/100)2 /(365 × 24)). Even if
we use 1000 nodes of the ABCI (which has 1088 nodes[51]) in parallel, SGLT(T) would require more than four
years. Consequently, SGLT cannot be applied to IST.
Therefore, we compare PPMTF with SGD in FS.

Uniformk = 100

0.25

TM-EMD-Y

Fig. 9. Running time in PF. “T” and “S” in the parentheses represent the time to train a generative model (i.e., MTF parametersin PPMTF and semantic clusters in SGLT) and the time to generate 500 synthetic traces, respectively.

VF-TV

|𝒰|

TP-TV-Top50|𝓧|=400

0.6
0.5
0.4
0.3
0.2
0.1
0

TP-TV

SGLT (S)

102
101
100
10-1
10-2
10-3
10-4
10-5
10-6

100
200
300
400
500
600
700
800
900
1000

Running Time (hour)|𝓤|=100|𝒳|

SGLT (T)PPMTF (S)

100
200
300
400
500
600
700
800
900
1000

Running Time (hour)

PPMTF (T)
102
101
100
10-1
10-2
10-3
10-4
10-5
10-6

14

IST

NYC

TKY

0

Fig. 10. Utility of synthetic traces with (k, 1)-PD in FS. The left
graphs show the utility of PPMTF without the PD test, withk = 10, 100 or 200. Lower is better in all of the utility metrics.

is missing for many users in FS, and cannot be copied.
Note that SGD with ξ = 0 always passes the PD test
because it generates synthetic traces independently of
the input data record [9]. We evaluated all the utility
metrics for PPMTF and SGD.
Figure 10 shows the results. The left graphs show
PPMTF without the PD test, with k = 10, 100, or 200
in IST, NYC, and TKY (we confirmed that the results of
the other cities were similar to those of NYC and TKY).
The right graphs show PPMTF with k = 10 and SGD.
The left graphs show that all of the utility metrics are minimally affected by running the PD test withk = 10 in all of the cities. Similarly, all of the utility metrics are minimally affected in IST, even whenk = 200. We confirmed that about 70% of the synthetic
traces passed the PD test when k = 10, whereas only
about 20% of the synthetic traces passed the PD testwhen k = 200 (see Appendix E for details). Nevertheless, PPMTF significantly outperforms Uniform in IST.
This is because the number of users is very large in
IST (|U| = 219793). Consequently, even if the PD test
pass rate is low, many synthetic traces still pass the test

15

Privacy-Preserving Multiple Tensor Factorization for Synthesizing Large-Scale Location Traces with Cluster-Specific Features

Frequency

Frequency(i) Park

1
Testing
0.8
PPMTF
0.6
SGD
0.4
0.2
0 1 2 3 4 5 6 7 8 9 10 11 12 13-

Visit-fraction (/24)

Frequency

Frequency(iii) Bar1
Testing
0.8
PPMTF
0.6
SGD
0.4
0.2
0 1 2 3 4 5 6 7 8 9 10 11 12 13-

Visit-fraction (/24)(i) Cluster of great outdoors and malls(ii) Office1
Testing
0.8
PPMTF
0.6
SGD
0.4
0.2
0 1 2 3 4 5 6 7 8 9 10 11 12 13-

Visit-fraction (/24)1(iv) University

Testing
0.8
PPMTF
0.6
SGD
0.4
0.2
0 1 2 3 4 5 6 7 8 9 10 11 12 13-

Visit-fraction (/24)

Fig. 11. Distributions of visit-fractions in NYC. PPMTF provides(10, 1)-PD.

Great Outdoors (0)Great Outdoors (14)Great Outdoors (3)Mall (15)Mall (9)01.6Analysis on Cluster-Specific Features. First, weshow in Figure 11 the distributions of visit-fractions forfour POI categories in NYC (Testing represents the distribution of testing traces). The distribution of SGD
concentrates at the visit-fraction of 1/24 (i.e., 0 to
0.042). This is because SGD (ξ = 0) uses the transition matrix and visit-probability vector common to
all users, and synthesizes traces independently of input
users. Consequently, all users spend almost the sameamount of time on each POI category. On the other
hand, PPMTF models a histogram of visited locations
for each user via the visit-count tensor, and generates
traces based on the tensor. As a result, the distributionof PPMTF is similar to that of Testing, and reflects thefact that about 30 to 35% of users spend less than 1/24of their time at a park or bar, whereas about 80% ofusers spend more than 1/24 of their time at an office orSynthetic (10-12h)

Great Outdoors (14)Great Outdoors (14)Mall (9)Great Outdoors (3)University (5)Mall (15)Great Outdoors (0)Great Outdoors (0)Great Outdoors (3)Airport (2)80 0.01 0.02
0 0.02 0.04

4

Training (10-12h)

Factor matrix D

Synthetic (10-12h)

0.8
00

6

12

Hour

18

24

(ii) Cluster of universities
Training (10-12h)

Factor matrix B
University (5)University (16)Great Outdoors (0)University (56)
University (33)
0

18

36

Great Outdoors (0)University (16)
University (32)
University (5)Mall (12)0 0.025 0.05

Training (10-12h)

Factor matrix D

Synthetic (10-12h)
Great Outdoors (0)University (16)
University (5)Bridge (6)Neighborhood (10)00.01 0.02

Synthetic (10-12h)

2.6
1.3
00

and preserve various statistical features. Thus PPMTF
achieves high utility especially for a large-scale dataset.
The right graphs in Figure 10 show that for TP-TV
and TP-TV-Top50, PPMTF is roughly the same as SGD.
For TM-EMD-X and TM-EMD-Y, PPMTF outperforms
SGD, especially in JK, NYC, KL, and SP. This is because
many missing events exist in FS and the transitions in
the training traces are few in JK, NYC, KL, and SP (asdescribed in Section 4.1).
A crucial difference between PPMTF and SGD lies in
the fact that PPMTF models the cluster-specific mobility features (i.e., both (c) and (d)), whereas SGD (ξ = 0)
does not. This causes the results of VF-TV in Figure 10.
Specifically, for VF-TV, SGD performs almost the sameas Uniform, whereas PPMTF significantly outperforms
SGD. Below we perform more detailed analysis to show
how well PPMTF provides (c) and (d).

Training (10-12h)

Factor matrix B

6

12

Hour

18

24

Fig. 12. Two clusters in IST (21980 users for each cluster). HerePPMTF provides (10, 1)-PD. For B and training/synthetic traces,
we show the top 5 POIs (numbers in parentheses represent POIIDs), whose values or frequencies from 10:00 to 12:00 are thehighest. We show the top 20 POIs by circles in the map. Red circles in (i) (resp. (ii)) represent outdoors/malls (resp. universities).

university. This result explains the low values of VF-TV
in PPMTF. Figure 11 also shows that PPMTF is usefulfor semantic annotation of POIs [19, 70].
Next, we visualize in Figure 12 the columns of factor matrices B and D and training/synthetic traces fortwo clusters. As with PF, the training users in eachcluster exhibit a similar behavior; e.g., the users in (i)enjoy great outdoors and shopping at a mall, whereas
the users in (ii) go to universities. Note that users and
POIs in each cluster are semantically similar; e.g., people who enjoy great outdoors also enjoy shopping at
a mall; many users in (ii) would be students, faculty,
or staff. The activity times are also different between
the two clusters. For example, we confirmed that many
training users in (i) enjoy great outdoors and shoppingfrom morning until night, whereas most training usersin (ii) are not at universities at night. PPMTF models such a behavior via factor matrices, and synthesizes
traces preserving the behavior. We emphasize that thisfeature is useful for various analysis; e.g., modeling human location patterns, semantic annotation of POIs.
SGD (ξ = 0) and others [12, 13, 28] do not provide such cluster-specific features because they generate
traces only based on parameters common to all users.
Scalability. Figure 13 shows the running time in FS.
SGD is much faster than PPMTF. The reason for thislies in the simplicity of SGD; i.e., SGD trains a transition

24
18
12
6
0

PPMTF (T)PPMTF (S)PPMTF (P)IST JK NYC KL SP TKY

100

Time (sec)Time (hour)

Privacy-Preserving Multiple Tensor Factorization for Synthesizing Large-Scale Location Traces with Cluster-Specific Features

75
50

SGD (T)
SGD (S)

25
0References[1]

IST JK NYC KL SP TKY

Fig. 13. Running time in FS. “T”, “S”, and “P” represent thetime to train a generative model, synthesize traces, and run the
PD test, respectively.

matrix for each time slot via maximum likelihood estimation; it then synthesizes traces using the transition
matrix. However, SGD does not generate cluster-specifictraces. To generate such traces, PPMTF is necessary.
Note that even though we used a supercomputerin our experiments, we used a single node and did notparallelize the process. We can also run PPMTF on aregular computer with large memory. For example, assume that we use 8 bytes to store a real number, andthat we want to synthesize all of 219793 traces in IST.
Then, 8|U|(λI + ρI + λII + ρII ) + 8z(|U| + 2|X | + |L|) = 3.9
GB memory is required to perform MTF, and the otherprocesses need less memory. PPMTF could also be parallelized by using asynchronous Gibbs sampling [66].

5 Conclusion
In this paper, we proposed PPMTF (Privacy-Preserving
Multiple Tensor Factorization), a location synthesizer
that preserves various statistical features, protects user
privacy, and synthesizes large-scale location traces inpractical time. Our experimental results showed that
PPMTF significantly outperforms two state-of-the-art
location synthesizers [8, 9] in terms of utility and scalability at the same level of privacy.
We assumed a scenario where parameters of the generative model are kept secret (or discarded after synthesizing traces). As future work, we would like to design alocation synthesizer that provides strong privacy guarantees in a scenario where the parameters of the generative model are made public. For example, one possibility
might be to release only parameters (B, C, D) (i.e., location and time profiles) and randomly generate A (i.e.,user profile) from some distribution. We would like toinvestigate how much this approach can reduce ε in DP.
Acknowledgement: This study was supported by JSPS
KAKENHI JP19H04113, JP17K12667, and by Inria under the project LOGIS.

16

Tool: Privacy-preserving multiple tensor factorization(PPMTF). https://github.com/PPMTF/PPMTF.
[2] PWS Cup 2019.
https://www.iwsec.org/pws/2019/cup19_e.html,
2019.
[3] C. C. Aggarwal. Recommender Systems. Springer, 2016.
[4] R. Albright, J. Cox, D. Duling, A. N. Langville, and C. D.
Meyer. Algorithms, initializations, and convergence for thenonnegative matrix factorization. SAS Technical Report,
pages 1–18, 2014.
[5] J. Biagioni and J. Eriksson. Inferring road maps from global
positioning system traces: Survey and comparative evaluation. Journal of the Transportation Research Board,2291(2291):61–71, 2012.
[6] I. Bilogrevic, K. Huguenin, M. Jadliwala, F. Lopez, J.-P.
Hubaux, P. Ginzboorg, and V. Niemi. Inferring social ties inacademic networks using short-range wireless communications. In Proc. WPES’13, pages 179–188, 2013.
[7] V. Bindschaedler and R. Shokri. Synthetic location traces
generator (sglt). https://vbinds.ch/node/70.
[8] V. Bindschaedler and R. Shokri. Synthesizing plausible
privacy-preserving location traces. In Proc. S&P’16, pages
546–563, 2016.
[9] V. Bindschaedler, R. Shokri, and C. A. Gunter. Plausible
deniability for privacy-preserving data synthesis. PVLDB,
10(5):481–492, 2017.
[10] C. Bishop. Pattern Recognition and Machine Learning.
Springer, 2006.
[11] K. Chatzikokolakis, E. Elsalamouny, C. Palamidessi, and
A. Pazii. Methods for location privacy: A comparative
overview. Foundations and Trends in Privacy and Security,
1(4):199–257, 2017.
[12] R. Chen, G. Acs, and C. Castelluccia. Differentially privatesequential data publication via variable-length n-grams. In
Proc. CCS’12, pages 638–649, 2012.
[13] R. Chen, B. C. M. Fung, B. C. Desai, and N. M. Sossou.
Differentially private transit data publication: A case study
on the montreal transportation system. In Proc. KDD’12,
pages 213–221, 2012.
[14] E. Cho, S. A. Myers, and J. Leskovec. Friendship and mobility: User movement in location-based social networks. In
Proc. KDD’11, pages 1082–1090, 2011.
[15] R. Chow and P. Golle. Faking contextual data for fun, profit,and privacy. In Proc. WPES’09, pages 105–108, 2009.
[16] A. Cichocki, R. Zdunek, A. H. Phan, and S. Amari. Nonnegative Matrix and Tensor Factorizations: Applications to
Exploratory Multi-way Data Analysis and Blind Source Separation. Wiley, 2009.
[17] J. Cranshaw, R. Schwartz, J. I. Hong, and N. Sadeh. The
livehoods project: Utilizing social media to understand the
dynamics of a city. In Proc. ICWSM’12, pages 58–65, 2012.
[18] Y.-A. de Montjoye, C. A. Hidalgo, M. Verleysen, and V. D.
Blondel. Unique in the crowd: The privacy bounds of humanmobility. Scientific Reports, 3(1376):1–5, 2013.
[19] T. M. T. Do and D. Gatica-Perez. The places of our lives:
Visiting patterns and automatic labeling from longitudinalsmartphone data. IEEE Trans. Mob. Comput, 13(3):638–

Privacy-Preserving Multiple Tensor Factorization for Synthesizing Large-Scale Location Traces with Cluster-Specific Features

648, 2013.
[20] J. Domingo-Ferrer, S. Ricci, and J. Soria-Comas. Disclosurerisk assessment via record linkage by a maximum-knowledge
attacker. In Proc. PST’15, pages 3469–3478, 2015.
[21] C. Dwork. Differential privacy. In Proc. ICALP’06, pages
1–12, 2006.
[22] C. Dwork and A. Roth. The Algorithmic Foundations ofDifferential Privacy. Now Publishers, 2014.
[23] C. Dwork and A. Smith. Differential privacy for statistics:
What we know and what we want to learn. Journal of Privacy and Confidentiality, 1(2):135–154, 2009.
[24] N. Eagle, A. Pentland, and D. Lazer. Inferring friendshipnetwork structure by using mobile phone data. PNAS,
106(36):15274–15278, 2009.
[25] J. Ernvall and O. Nevalainen. An algorithm for unbiased
random sampling. The Computer Journal, 25(1):45–47,
1982.
[26] S. Gambs, M.-O. Killijian, and M. Núñez del Prado Cortez.
De-anonymization attack on geolocated data. Journal ofComputer and System Sciences, 80(8):1597–1614, 2014.
[27] G. Ghinita. Privacy for Location-based Services. Morgan &Claypool Publishers, 2013.
[28] X. He, G. Cormode, A. Machanavajjhala, C. M. Procopiuc,and D. Srivastava. DPT: Differentially private trajectory
synthesis using hierarchical reference systems. PVLDB,
11(8):1154–1165, 2015.
[29] H. Hu, J. Xu, Q. Chen, and Z. Yang. Authenticating
location-based services without compromising location privacy. In Proc. SIGMOD’12, pages 301–312, 2012.
[30] T. Iwata and H. Shimizu. Neural collective graphical models
for estimating spatio-temporal population flow from aggregated data. In Proc. AAAI’19, pages 3935–3942, 2019.
[31] B. Jayaraman and D. Evans. Evaluating differentially private
machine learning in practice. In Proc. USENIX Security’19,
pages 1895–1912, 2019.
[32] R. Kato, M. Iwata, T. Hara, A. Suzuki, X. Xie, Y. Arase,
and S. Nishio. A dummy-based anonymization method
based on user trajectory with pauses. In Proc. SIGSPATIAL’12, pages 249–258, 2012.
[33] Y. Kawamoto and T. Murakami. Local obfuscation mechanisms for hiding probability distributions. In Proc. ESORICS,
pages 128–148, 2019.
[34] R. H. Keshavan, A. Montanari, and S. Oh. Matrix completion from noisy entries. In Proc. NIPS’09, pages 952–960,
2009.
[35] S. A. Khan and S. Kaski. Bayesian multi-view tensor factorization. In Proc. ECML PKDD’14, pages 656–671, 2014.
[36] H. Kido, Y. Yanagisawa, and T. Satoh. An anonymous
communication technique using dummies for location-based
services. Proc. ICPS’05, pages 88–97, 2005.
[37] J. Krumm. A survey of computational location privacy.
Personal and Ubiquitous Computing, 13(6):391–399, 2009.
[38] N. Li, M. Lyu, and D. Su. Differential Privacy: From Theory
to Practice. Morgan & Claypool Publishers, 2016.
[39] L. Liao, D. Fox, and H. Kautz. Extracting places and activities from gps traces using hierarchical conditional random fields. International Journal of Robotics Research,
26(1):119–134, 2007.
[40] M. Lichman and P. Smyth. Modeling human location datawith mixtures of kernel densities. In Proc. KDD’14, pages

17

35–44, 2014.
[41] X. Liu, J. Biagioni, J. Eriksson, Y. Wang, G. Forman, andY. Zhu. Mining large-scale, sparse gps traces for map inference: Comparison of approaches. In Proc. KDD’12, pages
669–677, 2012.
[42] X. Liu, Y. Liu, K. Aberer, and C. Miao. Personalized pointof-interest recommendation by mining users’ preference transition. In Proc. CIKM’13, pages 733–738, 2013.
[43] Z. Liu, Y.-X. Wang, and A. J. Smola. Fast differentiallyprivate matrix factorization. In Proc. RecSys’15, pages 171–
178, 2015.
[44] Y. Matsuo, N. Okazaki, K. Izumi, Y. Nakamura,T. Nishimura, and K. Hasida. Inferring long-term user properties based on users’ location history. In Proc. IJCAI’07,
pages 2159–2165, 2007.
[45] X. Meng, S. Wang, K. Shu, J. Li, B. Chen, H. Liu, andY. Zhang. Personalized privacy-preserving social recommendation. In Proc. AAAI’18, pages 1–8, 2018.
[46] T. Murakami. Expectation-maximization tensor factorization
for practical location privacy attacks. PoPETs, 4:138–155,
2017.
[47] T. Murakami, A. Kanemura, and H. Hino. Group sparsity
tensor factorization for de-anonymization of mobility traces.
In Proc. TrustCom’15, pages 621–629, 2015.
[48] T. Murakami, A. Kanemura, and H. Hino. Group sparsity
tensor factorization for re-identification of open mobility
traces. IEEE Trans. Inf. Forensics Secur., 12(3):689–704,
2017.
[49] T. Murakami and H. Watanabe. Localization attacks using
matrix and tensor factorization. IEEE Trans. Inf. Forensics
Secur., 11(8):1647–1660, 2016.
[50] K. P. Murphy. Machine Learning: A Probabilistic Perspective. The MIT Press, 2012.
[51] National Institute of Advanced Industrial Science and Technology (AIST). AI bridging cloud infrastructure (ABCI).
https://abci.ai/.
[52] Nightley and Center for Spatial Information Science at the
University of Tokyo (CSIS). SNS-based people flow data.
http://nightley.jp/archives/1954, 2014.
[53] V. Nikolaenko, S. Ioannidis, U. Weinsberg, M. Joye, N. Taft,and D. Boneh. Privacy-preserving matrix factorization. In
Proc. CCS’13, pages 801–812, 2013.
[54] R. Pan, Y. Zhou, B. Cao, N. N. Liu, R. Lukose, M. Scholz,
and Q. Yang. One-class collaborative filtering. In Proc.
ICDM’08, pages 502–511, 2008.
[55] M. Piorkowski, N. Sarafijanovic-Djukic, and M. Grossglauser. CRAWDAD dataset epfl/mobility (v. 2009-02-24).
http://crawdad.org/epfl/mobility/20090224, 2009.
[56] V. Primault, A. Boutet, S. B. Mokhtar, and L. Brunie. Thelong road to computational location privacy: A survey. IEEE
Commun. Surv., 21(3):2772–2793, 2019.
[57] R. Salakhutdinov and A. Mnih. Probabilistic matrix factorization. In Proc. NIPS’07, pages 1257–1264, 2007.
[58] R. Salakhutdinov and A. Mnih. Bayesian probabilistic matrix
factorization using markov chain monte carlo. In Proc.
ICML’08, pages 880–887, 2008.
[59] Y. Sekimoto, R. Shibasaki, H. Kanasugi, T. Usui, andY. Shimazaki. PFlow: Reconstructing people flow recycling
large-scale social survey data. IEEE Pervasive Computing,10(4):27–35, 2011.

Privacy-Preserving Multiple Tensor Factorization for Synthesizing Large-Scale Location Traces with Cluster-Specific Features[60] S. Shekhar, M. R. Evans, V. Gunturi, and K. Yang. Spatial
big-data challenges intersecting mobility and cloud computing. In Proc. MobiDE’12, pages 1–12, 2012.
[61] R. Shokri, M. Stronati, C. Song, and V. Shmatikov. Membership inference attacks against machine learning models.
In Proc. S&P’17, pages 3–18, 2017.
[62] R. Shokri, G. Theodorakopoulos, J.-Y. L. Boudec, and J.-P.
Hubaux. Quantifying location privacy. In Proc. S&P’11,
pages 247–262, 2011.
[63] L. Song, D. Kotz, R. Jain, and X. He. Evaluating next-cell
predictors with extensive wi-fi mobility data. IEEE Trans.
Mob. Comput, 5(12):1633–1649, 2006.
[64] A. Suzuki, M. Iwata, Y. Arase, T. Hara, X. Xie, andS. Nishio. A user location anonymization method for location based services in a real environment. In Proc. GIS’10,
pages 398–401, 2010.
[65] K. Takeuchi, R. Tomioka, K. Ishiguro, A. Kimura, and
H. Sawada. Non-negative multiple tensor factorization.
In Proc. ICDM’13, pages 1199–1204, 2013.
[66] A. Terenin, D. Simpson, and D. Draper. Asynchronous gibbs
sampling. In Proc. AISTATS’20, pages 144–154, 2020.
[67] Y.-X. Wang, S. E. Fienberg, and A. J. Smola. Privacy for
free: Posterior sampling and stochastic gradient monte carlo.
In Proc. ICML’15, pages 2493–2502, 2015.
[68] D. Yang, B. Qu, J. Yang, and P. Cudre-Mauroux. Revisiting user mobility and social relationships in LBSNs: A hypergraph embedding approach. In Proc. WWW’19, pages
2147–2157, 2019.
[69] D. Yang, D. Zhang, and B. Qu. Participatory cultural mapping based on collective behavior data in location based social network. ACM Trans. Intell. Syst. Technol., 7(3):30:1–
30:23, 2016.
[70] M. Ye, D. Shou, W.-C. Lee, P. Yin, and K. Janowicz. On
the semantic annotation of places in location-based social
networks. In Proc. KDD’11, pages 520–528, 2011.
[71] S. Yeom, I. Giacomelli, M. Fredrikson, and S. Jha. Privacy
risk in machine learning: Analyzing the connection to overfitting. In Proc. CSF’18, pages 268–282, 2018.
[72] T.-H. You, W.-C. Peng, and W.-C. Lee. Protecting moving
trajectories with dummies. In Proc. MDM’07, pages 278–
282, 2007.
[73] V. W. Zheng, Y. Zheng, and Q. Yang. Joint learning user’s
activities and profiles from GPS data. In Proc. LBSN’09,
pages 17–20, 2009.
[74] Y. Zheng, X. Xie, and W.-Y. Ma. GeoLife: A collaborativesocial networking service among user, location and trajectory. IEEE Data Engineering Bulletin, 32(2):32–40, 2010.
[75] Y. Zheng, L. Zhang, X. Xie, and W.-Y. Ma. Mining interesting locations and travel sequences from GPS trajectories. In
Proc. WWW’09, pages 791–800, 2009.

A Notations and Abbreviations
Tables 1 and 2 respectively show the basic notations andabbreviations used in this paper.

18

Table 1. Basic notations in this paper († represents I or II).
Symbol
U
X
T
L
E
R
S
F
M
un
xi
sn
y
R
R̂†
†
rn,i,jDescription
Finite set of training users.
Finite set of locations.
Finite set of time instants over N.
Finite set of time slots (L ⊆ P(T )).
Finite set of events (E = X × T ).
Finite set of traces (R = U × E ∗ ).
Finite set of training traces (S ⊆ R).
Randomized algorithm with domain P(R).
Generative model.
n-th training user (un ∈ U ).
i-th location (xi ∈ X ).
n-th training trace (sn ∈ S).
Synthetic trace (y ∈ R).
Tuple of two tensors (R = (RI , RII )).
Reconstructed tensors by Θ.
(n, i, j)-th element of R† .

r̂n,i,jΘz
FPPMTFMPPMTFQn,i(n, i, j)-th element of R̂† .
Tuple of MTF parameters (Θ = (A, B, C, D)).
Number of columns in each factor matrix.
Proposed training algorithm.
Proposed generative model.
Transition-probability matrix of user un fortime slot li in MPPMTF .
Visit-probability vector of user un for timeslot li in MPPMTF .
Maximum number of positive elements peruser in R† .
Number of selected zero elements per userin R† .
Maximum value of counts for each elementin R† .
†
Indicator function that takes 0 if rn,i,j ismissing, and takes 1 otherwise.

†

πn,i
λ†
ρ†
†

rmax
†

In,i,jB Time Complexity
Assume that we generate a synthetic trace from eachtraining trace sn ∈ S (i.e., |U| synthetic traces in total).
Assume that λI , ρI , λII , ρII , z, and |U ∗ | are constants.
In step (i), we simply count the number of transitions and the number of visits from a training trace set
S. Consequently, the computation time of this step ismuch smaller than that of the remaining three steps.
In step (ii), we first randomly select ρI and ρII zeroelements for each user in RI and RII , respectively. Thiscan be done in O(|U|) time in total by using a samplingtechnique in [25]. Subsequently, we train the MTF parameters Θ via Gibbs sampling. The computation timeof Gibbs sampling can be expressed as O(|U|+|X |+|L|).
In step (iii), we generate synthetic traces via theMH algorithm. This is dominated by computation of thetransition-probability matrices Q∗n , Qn,1 , · · · , Qn,|L| forPrivacy-Preserving Multiple Tensor Factorization for Synthesizing Large-Scale Location Traces with Cluster-Specific Features

Table 2. Abbreviations in this paper.
Abbreviation
PPMTF
SGLT
SGD
PF
FS
IST/JK/NYC/
KL/SP/TKY
TP-TV(-Top50)

TM-EMD-X/Y
VF-TV

Description
Proposed location traces generator.
Synthetic location traces generator in [8].
Synthetic data generator in [9].
SNS-based people flow data [52].
Foursquare dataset [68].
Istanbul/Jakarta/New York City/
Kuala Lumpur/San Paulo/Tokyo.
Average total variation between timedependent population distributions(over 50 frequently visited locations).
Earth Mover’s Distance between transitionprobability matrices over the x/y-axis.
Total variation between distributions ofvisit-fractions.

each training trace sn , which takes O(|U||X |2 |L|) timein total. Then we generate a synthetic trace y, whichtakes O(|U||X ||L|) time.
In step (iv), the faster version of Privacy Test 1 in
Section 3.5 computes the transition-probability matrices
Q∗m , Qm,1 , · · · , Qm,|L| for each training trace sm ∈ S ∗ ,which takes O(|X |2 |L|) time in total. Subsequently, we
check whether k 0 ≥ k for each training trace sn ∈ S,
which takes O(|U||X ||L|) time in total.
In summary, the time complexity of the proposed
method can be expressed as O(|U||X |2 |L|).

C Details on SGD
SGD [9] is a synthetic generator for any kind of data,
which works as follows: (i) Train the dependency structure (graph) between data attributes; (ii) Train conditional probabilities for each attribute given its parent
attributes; (iii) Generate a synthetic data record froman input data record by copying the top γ ∈ Z≥0 attributes from the input data record and generating theremaining attributes using the trained conditional probabilities. Note that the dependency structure and theconditional probabilities are common to all users.
We applied SGD to synthesis of location traces asfollows. We regarded an event as an attribute, and alocation trace of length |T | as a data record with |T | attributes. Then it would be natural to consider that thedependency structure is given by the time-dependent
Markov chain model as in PPMTF and SGLT, and theconditional probabilities are given by the transition matrix for each time slot. In other words, we need not train
the dependency structure; i.e., we can skip (i).

19

We trained the transition matrix Q̃i ∈ Q for eachtime slot li ∈ L (|L| × |X | × |X | elements in total) andthe visit-probability vector π̃ ∈ C for the first time instant (|X | elements in total) from the training traces
via maximum likelihood estimation. Then we synthesized a trace from an input user un by copying the firstγ events in the training trace sn of un and by generating
the remaining events using the transition matrix. Whenγ = 0, we generated a location at the first time instantusing the visit-probability vector. Thus the parametersof the generative model Mn of user un can be expressed
as: (Q̃1 , · · · , Q̃|L| , π̃, sn ).
SGD can provide (ε, δ)-DP for one synthetic trace y
by using a randomized test [9], which randomly selects
an input user un from U and adds the Laplacian noise
to the parameter k in (k, η)-PD. However, both ε andδ can be large for multiple synthetic traces generated
from the same input user, as discussed in [9]. Thus we
did not use the randomized test in our experiments.

D Details on Privacy Attacks
Re-Identification Algorithm. We used the Bayesian
re-identification algorithm in [47]. Specifically, we firsttrained the transition matrix for each training user fromthe training traces via maximum likelihood estimation.
Then we re-identified each synthetic trace y by selecting
a training user whose posterior probability of being the
input user is the highest. Here we computed the posterior probability by calculating a likelihood for eachtraining user and assuming a uniform prior for users.
We calculated the likelihood by simply calculating alikelihood for each transition in y using the transition
matrix and multiplying them. We assigned a small positive value (= 10−8 ) to zero elements in the transition
matrix so that the likelihood never becomes 0.
Membership Inference Algorithm. We considered a likelihood ratio-based membership inference algorithm, which partly uses the algorithm in [48] as follows.
Let V be a finite set of all training and testing users(each of them is either a member or a non-member;|V| = 1000 in PF), and vn ∈ V be the n-th user. Assume that the adversary attempts to determine whether
user vn is a training user (i.e., member) or not. Since
each training user is used as an input user to generate asynthetic trace, the adversary can perform the membership inference by determining, for each synthetic tracey, whether vn is used as an input user to generate y.
To perform this two-class classification (i.e., vn is an in-put user of y or not), we used the likelihood ratio-based
two-class classification algorithm in [48].
Specifically, given user vn and synthetic trace y, let
H1 (resp. H0 ) be the hypothesis that vn is (resp. is not)an input user of y. We first trained the transition matrix
for each of |V| = 1000 users from her (original) trace.
Let Wn be the transition matrix of user vn . We calculated a population transition matrix W0 , which models
the average behavior of users other than vn as the avP
1
erage of Wm (m 6= n); i.e., W0 = |V|−1
m6=n Wm .
Let z1 (resp. z0 ) ∈ R be the likelihood of y given H1(resp. H0 ). We calculated z1 (resp. z0 ) simply by calculating a likelihood for each transition in y using thetransition matrix Wn (resp. W0 ) and multiplying them(as in the re-identification attack). Then we compared
the log-likelihood ratio log zz01 with a threshold ψ ∈ R. If
log zz10 ≥ ψ, we accepted H1 ; otherwise, we accepted H0 .
We performed this two-class classification for eachsynthetic trace y. If we accepted H1 for at least one
synthetic trace y, then we decided that vn is a member.
Otherwise, we decided that vn is a non-member.
In our experiments, we changed the threshold ψ tovarious values. Then we evaluated, for each location synthesizer, the maximum membership advantage over thevarious thresholds (Figure 8 shows the results).

E Relationship between k and the
PD Test Pass Rate
We evaluated the PD test pass rate, which is the proportion of synthetic traces that have passed the PD test
to all synthetic traces when we changed k from 1 to 200.
We set the other parameters to the same values as in
Section 4 (e.g., η = 1, |U ∗ | = 32000).
Figure 14 shows the results obtained for six cities
in FS. The PD test pass rate decreases with an increasein k. For example, the PD test pass rate is about 70%
when k = 10, whereas it is about 20% when k = 200.
Note that when k = 200, the PD test pass rateof IST (17.9%) is lower than that of NYC (26.9%), as
shown in Figure 14. Nevertheless, PPMTF significantly
outperforms Uniform with regard to all of the utility
metrics in IST, as shown in Figure 10. This is because
the number of users is very large in IST (|U| = 219793).
Consequently, even if the PD test pass rate is low, many
synthetic traces still pass the test and preserve various
statistical features.
Therefore, PPMTF achieves high utility especiallyfor a large-scale dataset.

PD Test Pass Rate

Privacy-Preserving Multiple Tensor Factorization for Synthesizing Large-Scale Location Traces with Cluster-Specific Features

1

IST
NYC
SP

0.8
0.6

20

JK
KL
TKY

0.4
0.2
0

0

50

100
k

150

200

Fig. 14. Relationship between k and the PD test pass rate.

F DP for the MTF Parameters ΘHere we explain DP (Differential Privacy) [21, 22] asa privacy metric (Appendix F.1). Then we analyze theprivacy budget ε in DP for the MTF parameters Θ in
PPMTF (Appendix F.2), and evaluate ε for Θ using theFoursquare dataset (Appendix F.3).

F.1 Differential Privacy
We define the notion of neighboring data sets in thesame way as [22, 43, 67] as follows. Let S, S 0 ⊆ R be two
sets of training traces. We say S and S 0 are neighboring
if they differ by at most one trace and include the same
number of traces, i.e., |S| = |S 0 |. For example, given a
trace s01 ∈ R, S = {s1 , s2 , s3 } and S 0 = {s01 , s2 , s3 } areneighboring. Then DP [21, 22] is defined as follows:
Definition 2 (ε-DP). Let ε ∈ R≥0 . A randomized algorithm F with domain P(R) provides ε-DP if for any
neighboring S, S 0 ⊆ R and any Z ⊆ Range(F),e−ε p(F(S 0 ) ∈ Z) ≤ p(F(S) ∈ Z) ≤ eε p(F(S 0 ) ∈ Z). (8)ε-DP guarantees that an adversary who has observed
the output of F cannot determine, for any pair of S
and S 0 , whether it comes from S or S 0 (i.e., a particularuser’s trace is included in the training trace set) witha certain degree of confidence. As the privacy budget
ε approaches 0, S and S 0 become almost equally likely,
which means that a user’s privacy is strongly protected.

F.2 Theoretical Analysis
We now analyze the privacy budget ε in DP for the
MTF parameters Θ in PPMTF.
Let FPPMTF be our training algorithm in step (ii),which takes as input the training trace set S and outputs
the MTF parameters Θ. Assume that Θ is sampled fromthe exact posterior distribution p(Θ|R).

Privacy-Preserving Multiple Tensor Factorization for Synthesizing Large-Scale Location Traces with Cluster-Specific Features

Recall that the maximum counts in RI and RII are
I
II , respectively, as defined in Section 3.2.
rmax
and rmax
Let κ ∈ R≥0 be a non-negative real number such that
I
I
II
II
r̂n,i,j∈ [−κ, rmax
+ κ] and r̂n,i,j∈ [−κ, rmax
+ κ] foreach triple (n, i, j). The value of κ can be made small
by iterating the sampling of Θ until we find Θ with smallκ [43]. Note that this “retry if fail” procedure guarantees
that Θ is sampled from the posterior distribution underIIIthe constraint that r̂n,i,j
and r̂n,i,j
are bounded as above(see the proof of Theorem 1 in [43]). Then we obtain:
Proposition 1. FPPMTF provides ε-DP, whereIε =α min{3λI , λI + ρI }(rmax
+ κ)2


II
+ min{3λII , λII + ρII }(rmax+ κ)2 .

(9)

Proof. By (3), ln p(Θ|R) can be written as follows:ln p(Θ|R)= ln p(R|Θ) + ln p(Θ) − ln p(R)(by Bayes’ theorem)r !
IIα(rn,i,j− r̂n,i,j
)2
α
+ ln
2
2π
r !
II
II )2
α(rn,i,j− r̂n,i,j
α
+ ln
2
2π

|U | |X | |X |

=−

XXX

I
In,i,jn=1 i=1 j=1−|U | |X | |L|
X
XX

II
In,i,jn=1 i=1 j=1+ ln p(Θ) − ln p(R).

(10)

The sum of the first and second terms in (10) is the loglikelihood ln p(R|Θ), and is bounded by the trimming
I
I
II
II ].
that ensures rn,i,j∈ [0, rmax
] and rn,i,j∈ [0, rmax
Let G be a function that takes as input R and Θ
and outputs G(R, Θ) ∈ R as follows:
G(R, Θ) =|U | |X | |X |IIIX
X X α In,i,j(rn,i,j− r̂n,i,j)2
n=1 i=1 j=1+

2|U | |X | |L|
II (r IIII
2X
X X α In,i,j
n,i,j − r̂n,i,j )n=1 i=1 j=1

2− ln p(Θ).
(11)
pαNote that ln 2π and ln p(R) in (10) do not depend on
Θ. Thus, by (11), ln p(Θ|R) in (10) can be expressed as:
exp[−G(R, Θ)].
exp[−G(R, Θ)]dΘΘ

p(Θ|R) = R(12)Then, Proposition 1 can be proven by using the factthat FP P M T F is the exponential mechanism [22] thatuses −G(R, Θ) as a utility function. Specifically, let R0
be the tuple of two tensors that differ from R at mostone user’s elements; i.e., R and R0 are neighboring. We

21write R ∼ R0 to represent that R and R0 are neighboring. Let ∆G ∈ R be the sensitivity of G given by:∆G = maxmaxΘ R,R0 :R∼R0|G(R, Θ) − G(R0 , Θ)|.

(13)

Here we note that when ρI is large, many zero elementsI
are common in RI and R0 . Specifically, for each user,
we can randomly select ρI zero elements as follows: (i)randomly select ρI elements from RI (including non-zeroelements), (ii) count the number ρI0 (≤ λI ) of non-zeroelements in the selected elements, (iii) randomly reselectρI0 elements from zero (and not selected) elements in
RI . Note that this algorithm eventually selects ρI zeroelements from RI at random.1 In this case, for each user,
at least max{ρI −2λI , 0} zero elements are common in RII
I
and R0 (since RI and R0 have at most 2λI reselected
elements in total).
I
Except for such common zero elements, In,i,jinIII(11) takes 1 at most min{3λ , λ + ρ } elements foreach user (since (λI + ρI ) − (ρI − 2λI ) = 3λI ). SimiII
larly, except for common zero elements, In,i,jin (11)II
II
II
takes 1 at most min{3λ , λ + ρ } elements for each
I
I
II
II ],
user. In addition, rn,i,j∈ [0, rmax], rn,i,j∈ [0, rmax
I
I
II
II +κ] for eachr̂n,i,j∈ [−κ, rmax+κ], and r̂n,i,j∈ [−κ, rmax
triple (n, i, j), as described in Section 3.5. Moreover, the“retry if fail” procedure, which iterates the sampling ofIIIΘ until r̂n,i,j
and r̂n,i,j
are bounded as above, guarantees that Θ is sampled from the posterior distribution
under this constraint [43].
Consequently, the sum of the first and second termsin (11) is less than (resp. more than) or equal to 2ε(resp. 0), where ε is given by (9). Then, since the third
term in (11) is the same for G(R, Θ) and G(R0 , Θ) in(13), ∆G can be bounded above by 2ε : i.e., ∆G ≤ 2ε .
Since the exponential mechanism with sensitivity 2ε provides ε-DP [22], FP P M T F provides ε-DP.
ε for a Single Location. We also analyze ε for neighboring data sets R and R0 that differ in a single location.
Here we assume ρI = ρII = 0 to simplify the analysis (ifρI > 0 or ρII > 0, then ε will be larger because selected
zero elements can be different in R and R0 ). In thisIIIcase, RI and R0 (resp. RII and R0 ) differ in at mosttwo (resp. four) elements, and the value in each element
differs by 1.2 Then by (11) and (13), we obtain:αII∆G ≤ 2 ·(rmax+ κ)2 − (rmax+ κ − 1)2
2

1 Other random sampling algorithms do not change our conclusion because p(Θ|S) is obtained by marginalizing R = (RI , RII ).
2 In RII and R0 II , we can consider the case where one
transition-count differs by 2 and two transition-counts differ by

TP-TV-Top50

2.5
2
1.5
1
0.5
0 -2
10

PPMTF
Uniform
Training
1102

�

PPMTF (X)
Uniform (X)Training (X)

1

102

�

104

106

PPMTF (Y)Uniform (Y)Training (Y)

104

106

108

VF-TV

TP-TV

0.6
0.5
0.4
0.3
0.2
0.1
0 -2
10

TM-EMD-X/YPrivacy-Preserving Multiple Tensor Factorization for Synthesizing Large-Scale Location Traces with Cluster-Specific Features

108

0.25
0.2
0.15
0.1
0.05
0 -2
10
0.6
0.5
0.4
0.3
0.2
0.1
0 -2
10

PPMTF
Uniform
Training
1

102

�

104

106

108

PPMTF
Uniform
Training

1

102

�

104

106

108

Fig. 15. Relation between ε and utility in IST (κ = 0).


α
II
II
(rmax
+ κ)2 − (rmax+ κ − 1)2 )2III+ 6κ − 3),+ 4rmax= α(2rmax+4·III
and therefore ε = α(4rmax+ 8rmax+ 12κ − 6).

Note that a trace y is synthesized from Θ after
FPPMTF outputs Θ. Then by the immunity to postprocessing [22], FPPMTF also provides ε-DP for all synthetic traces. However, ε needs to be large to achieve
high utility, as shown in Appendix F.3.

F.3 Experimental Evaluation
We evaluated the privacy budget ε in DP for Θ and theutility by changing α in Proposition 1 from 10−6 and
103 using the Foursquare dataset [68]. Figure 15 shows
the results in IST (Istanbul), where ε is the value in
Proposition 1 when κ = 0. In practice, ε can be larger
than this value because κ ≥ 0.
Figure 15 shows that ε needs to be larger than
2 × 104 to provide high utility. This is because α in
Proposition 1 needs to be large to achieve high utility.
I
Specifically, by (3), α needs to be large so that r̂n,i,j
and
II
I
II
r̂n,i,j in (3) are close to rn,i,j and rn,i,j , respectively. For
example, when α = 0.01 (i.e., standard deviation in (3)= 10), transition/visit-counts can be frequently changedII
by ±10 after sampling (e.g., r̂n,i,j= rn,i,j± 10), which
destroys the utility. In Figure 15, we need α ≥ 0.4 toachieve high utility, which results in ε > 2 × 104 .
If we consider neighboring data sets S and S 0 thatdiffer in a single location (rather than one trace), ε becomes much smaller. However, ε is still large. SpecifiIII
cally, if κ = ρI = ρII = 0, then ε = α(4rmax+ 8rmax−6);otherwise, ε is larger than this value (see Appendix F.2).

1 (e.g., transition x1 → x1 → x1 changes to x1 → x2 → x1 ). We
can ignore such cases because |G(R, Θ) − G(R0 , Θ)| is smaller.

22

Thus, when α = 0.4, the privacy budget is ε = 45.6 orIIImore (since rmax
= rmax= 10).
Finally, we note that adding the Laplacian noiseto Θ (rather than sampling Θ) does not provide DP.
For example, assume that Θ is trained from S by the
MAP (Maximum a Posteriori) estimation algorithm F[10], which calculates Θ that maximizes p(Θ|S); i.e.,
F(S) = argmaxΘ p(Θ|S). If p(Θ|S) is uniform (or nearly
uniform), then F(S 0 ) can take any value for neighboring
trace set S 0 . Therefore, the sensitivity is unbounded andadding the Laplacian noise does not provide DP.
For these reasons, providing a small ε in DP is difficult in our location synthesizer.