CiNCT: Compression and retrieval for massive
vehicular trajectories via relative movement labeling
Satoshi Koide #∗1 , Yukihiro Tadokoro #2 , Chuan Xiao ∗3 , Yoshiharu Ishikawa ∗4
#

arXiv:1706.02885v2 [] 29 Sep 2017

1

Toyota Central R&D Labs., Inc., Nagakute, Aichi, Japan
∗
Nagoya University, Nagoya, Aichi, Japan

koide@mosk.tytlabs.co.jp 2 tadokoro@mosk.tytlabs.co.jp
3
chuanx@nagoya-u.jp 4 ishikawa@i.nagoya-u.ac.jp

Abstract—In this paper, we present a compressed data structure for moving object trajectories in a road network, which
are represented as sequences of road edges. Unlike existing
compression methods for trajectories in a network, our method
supports pattern matching and decompression from an arbitrary
position while retaining a high compressibility with theoretical
guarantees. Specifically, our method is based on FM-index, a fast
and compact data structure for pattern matching. To enhance the
compression, we incorporate the sparsity of road networks into
the data structure. In particular, we present the novel concepts of
relative movement labeling and PseudoRank, each contributing to
significant reductions in data size and query processing time.
Our theoretical analysis and experimental studies reveal the
advantages of our proposed method as compared to existing
trajectory compression methods and FM-index variants.

I. I NTRODUCTION
In recent years, a vast amount of trajectory data from
moving objects, such as automobiles, has become available.
According to Han et al. [1], the total amount of GPS trajectories generated by automobiles in the U.S. alone exceeded
53 TB in 2011. With recent increased interest in the use of such
large datasets in wide range of data-driven applications, fundamental data manipulations such as retrieval and compression
are once again becoming crucial. In this paper, we focus on
moving object trajectories in (road) networks, called networkconstrained trajectories (NCTs), one of the most important
types of trajectories with many practical applications. Traveled
paths of NCTs can be represented as symbol sequences of road
segment IDs. Although this representation is more compact
than GPS coordinates, it is still insufficient for the vast datasets
that are now available. Therefore, compressed representations
of NCTs have been studied thus far [1–5].
If trajectories are simply compressed without an augmented
data structure, it is difficult to use them in real applications.
Therefore, compression methods that allow several operations
without decompressing the entire dataset are necessary, and
such methods have been the focus of recent studies. For
example, such studies include the in-memory data structures
proposed in [3] and [4], as well as an in-memory/on-disk
hybrid structure proposed in [6]. In our present paper, we
propose a method that realizes a high level of compression
while retaining a high utility of the data. As motivation and
background for our method, we first review existing compressed data structures for NCTs and their functions below.

NCTs consist of spatial paths and corresponding timestamps. We therefore must consider compression of these
paths and timestamps separately. For spatial paths, lossless
compression methods based on shortest-path encoding have
been studied in [1], [2], and [4]. Here, to compress the data,
these methods remove partial shortest paths in an NCT because
these paths can be recovered from the road network itself.
One drawback of this approach is that it cannot guarantee
the information-theoretic upper bound of the compressed data
size. A recent lossless path compressor introduced in [1]
called minimum entropy labeling (MEL) guarantees a theoretic
bound and also achieves practically higher compressibility
than shortest-path encoding methods. As for the timestamps,
all methods noted above compress them independently from
the spatial path compression. In this paper, we do not discuss
the compression of timestamps directly, but we emphasize here
that our method can be easily combined with such temporal
compression methods (see Section VII for details).
In general, it is difficult to define high utility of compressed
NCTs, because their utility depends on the given application.
In this paper, we focus on two functions, i.e., pattern matching
without decompressing the entire dataset, and extracting subpaths from an arbitrary position. Intuitively, pattern matching
operations that find trajectories along a given path would have
wide applications in NCT processing. In fact, the existing
methods mentioned above (i.e., [3], [4], and [6]) closely relate
to pattern matching; however, to the best of our knowledge,
there are no NCT compressors that guarantee theoretical bound
for the compressed size while supporting fast pattern matching.
Given the above, our research question is, how can we
realize high compressibility while enabling pattern matching
for NCTs? To address this question, we focus on suffix arrays
[7], data structures that closely relate to pattern matching.
Although the data structures for NCTs proposed in [3] and [6]
also employ suffix arrays, they do not focus on a compression
method, instead using existing general-purpose compressed
suffix arrays that are typically used for handling genomic sequences. Unfortunately, these existing methods are inefficient
because genomic sequences include only four characters (i.e,
A, C, G, and T) whereas NCTs consist of a large alphabet
(i.e., road segment IDs in a potentially large road network).
NCTs have another noteworthy feature, i.e., they can only
move along physically connected road segments. This feature

is quite different from general sequences, as illustrated in
Fig. 1. In Fig. 1(a), we show four example NCTs in a small
network with six road segments (A–F). The corresponding
graph shown in Fig. 1(b) represents symbol transitions for
these four NCTs. Here, each vertex corresponds to a symbol
(i.e., a road segment), and directed edges exist between two
vertices if the corresponding two symbols can appear successively. For example, in Fig. 1(b), vertex A is connected with
vertexes B and D because we can only move to road segment
B or D from A. For NCTs, this empirical transition graph
(ET-graph) becomes a sparse graph, reflecting the physical
topology of road networks. This sparsity cannot be obtained
for general sequences, which leads to a denser ET-graph, as
illustrated in Fig. 1(c).
Our proposed method, Compressed-index for NCTs
(CiNCT ), significantly improves the compression and pattern
matching operations when applied to sequences with such
sparse ET-graphs. Our method is based on FM-index [8], a
compressed data structure for suffix arrays, which we describe
further in Section II. Note that it is challenging to incorporate
such sparsity into FM-index while retaining its theoretical
advantages because FM-index is compressed at the bit-level.
Therefore, in the remainder of our paper, we introduce some
novel techniques and provide theoretical analysis that explains
why our method yields substantial improvement in practice.
Contributions: To develop a data structure for NCTs that
simultaneously achieves a high compression ratio and high
utility, we propose CiNCT, as a novel method to compress
suffix arrays for sequences on a sparse graph. We summarize
our contributions as follows.
• We propose relative movement labeling (RML), which
converts sequences on a sparse graph to low-entropy
sequences. We theoretically prove its optimality and show
that RML provides a more compact representation of
NCTs than that of the MEL method [1].
• We incorporate RML into FM-index by introducing a new
concept called PseudoRank, which leads to significant improvements in both size and query processing speed (i.e.,
the speed of pattern matching and sub-path extraction) as
compared to existing FM-index variants. We also explain
theoretically why this occurs.
• Using several real NCT datasets, we show that our
method outperforms state-of-the-art methods that do not
consider graph sparsity.
Outline: The remainder of our paper is organized as follows: preliminaries (Section II), proposed data structure (Section III), proposed algorithms (Section IV), theoretical analysis
(Section V), experiments (Section VI), related work (Section VII), and conclusion (Section VIII).
II. P RELIMINARIES
In this section, we introduce the data models and pattern
matching query. For readers not familiar with string processing and indexing, we also describe the necessary concepts
regarding FM-index and its compression. Table I summarizes
notation used in this paper.

A

B

C

B
E

D

C

A

D

F

F

(a) Network-constrained
trajectories (NCTs)

E

(b) ET-graph for
NCTs

(c) ET-graph for
general strings

Fig. 1. (a) Network-constrained trajectories (NCTs), and both (b) sparse and
(c) dense symbol transition graphs (ET-graphs).
TABLE I
N OTATION
Symbol
w, w0 ∈ E
T, Tbwt
Σ, σ
R(P ) = [sp, ep)
C[w]
H0 (S), Hk (S)
GT , ET
φ
Zw0 w

Description
Road segments (characters)
Trajectory string and its BWT
Alphabet set and its size
Suffix range of a pattern P
The number of w0 in T s.t. w0 < w
0th and kth order empirical entropy
ET-graph and its edge set
Relative movement labeling func.
Correction term

Defined in
—
Def. 2, Fig. 2
§ II-A1
§ II-A2
§ II-A3
Eq. (3), (4)
§ III-B (Def. 3)
§ III-B1
Eq. (7)

A. Definitions
1) Data models: First, we define NCTs as follows.
Definition 1: A network-constrained trajectory (NCT) on a
directed graph (V, E) is defined as a sequence of physically
connected road segments, i.e., e1 e2 · · · en (ei ∈ E).
For example, we have e1 = A, e2 = B, e3 = E, and e4 = F
for T1 = ABEF illustrated in Fig. 1 (a). To build an FM-index
for a set of documents, they are usually concatenated into one
long string [6]. Similarly, we define a trajectory string that
concatenates the NCTs.
Definition 2 (Trajectory string): Let T := {Tk }N
k=1 be a
set of NCTs to be indexed. A trajectory string is defined as
T := T1r $T2r $ · · · TNr $#, where Tkr is the reversal of string Tk ,
and $ and # are special symbols that represent NCT boundaries
and the end of the string, respectively.
For the four NCTs in Fig. 1 (a), the trajectory string is
T = FEBA
CB $ |{z}
DA $#.
| {z } $ |CBA
{z } $ |{z}
T1r

T2r

T3r

(1)

T4r

In the later sections, we use this example for explanation. In
this paper, a string S has 0-based subscripts and |S| denotes its
length. S[i] and S[i, j) are the i-th element and the substring
from i to j − 1, respectively. The alphabet set is defined
as Σ := E ∪ {$, #}, and σ denotes its size. To define the
BWT below, we assume a lexicographical order on the road
segments E (any ordering can be used for our purpose). The
lexicographical order is assumed to be # < $ < w (∀w ∈ E).
2) Pattern matching and BWT: The Burrows–Wheeler
transform (BWT) [9] is closely related to pattern matching
and is used in FM-index. It is a reversible transform of T ,
defined to be the last column of the lexicographically sorted
rotations of T (Fig. 2). For trajectory string Eq. (1), we have
Tbwt = $AAABDBBCCE$$$F#.

(2)

For a given pattern (string) P , we can define a unique range

All rotations of T

BWT:

Sorted rotations

Original string

T = FEBA$CBA$CB$DA$#

0:
1:
EBA$CBA$CB$DA$#F
2:
BA$CBA$CB$DA$#FE
3:
A$CBA$CB$DA$#FEB
4:
$CBA$CB$DA$#FEBA
5:
CBA$CB$DA$#FEBA$
BA$CB$DA$#FEBA$C Sort 6:
7:
A$CB$DA$#FEBA$CB
8:
$CB$DA$#FEBA$CBA
9:
CB$DA$#FEBA$CBA$
B$DA$#FEBA$CBA$C
10:
$DA$#FEBA$CBA$CB
BA 11:
12:
DA$#FEBA$CBA$CB$
13:
A$#FEBA$CBA$CB$D
14:
$#FEBA$CBA$CB$DA
#FEBA$CBA$CB$DA$
15:

#FEBA$CBA$CB$DA$
$#FEBA$CBA$CB$DA
$CB$DA$#FEBA$CBA
$CBA$CB$DA$#FEBA
$DA$#FEBA$CBA$CB
A$#FEBA$CBA$CB$D
A$CB$DA$#FEBA$CB
A$CBA$CB$DA$#FEB
B$DA$#FEBA$CBA$C
BA$CB$DA$#FEBA$C
BA$CBA$CB$DA$#FE
CB$DA$#FEBA$CBA$
CBA$CB$DA$#FEBA$
DA$#FEBA$CBA$CB$
EBA$CBA$CB$DA$#F
FEBA$CBA$CB$DA$#

Last
Column

$
A
A
A
B
D
B
B
C
C
E
$
$
$
F
#

Algorithm 1: Finding the suffix range R(P ) = [sp, ep) for a
given query P of length m based on Tbwt (SearchFM)

1
2
3
4
5
6
7

Fig. 2. The BWT of T is defined to be the last column of the sorted rotations
of T . This example is based on the trajectory string T in Eq. (1).

R(P ) = [sp, ep) for which the prefixes of the corresponding
sorted rotations are equal to P . We call this range the suffix
range of P . For example, if P = BA, we have R(P ) =
[9, 11) (see the underlined prefixes in Lines 9 and 10 in Fig. 2).
Finding R(P ) for a given P is called pattern matching, or
suffix range query, in this paper. The suffix range query for a
trajectory string T finds a suffix range of a given spatial path
P . It is known that the suffix range is useful in spatio-temporal
query processing for NCTs (see Section VII). This is why we
focus on this query.
In this paper, we also focus on another query, sub-path
extraction query, that recovers a sub-path of any length from
an arbitrary position in BWT Tbwt . We describe this query in
Section IV-C.
3) FM-index and an algorithm to find suffix ranges: The
FM-index [8] is a data structure that compresses a large string
and indexes it at the same time. Specifically, the FM-index of
a string T is a data structure in which the BWT of T is stored
in a wavelet tree. Suffix range queries can be processed rapidly
by FM-index. In the following, we overview how it works. It is
known that Algorithm 1 can find the suffix range R(P ) for any
P based on Tbwt . The rank function, rankw (Tbwt , i), returns
the number of occurrences of a symbol w ∈ Σ in a substring
Tbwt [0, i). For example, we have rankB (Tbwt , 5) = 1 because
Tbwt [0,5)

Tbwt

z }| {
= $AAAB DBBCCE$$$F#.

Moreover, C[w] is the number of symbols in Tbwt that
are lexicographically smaller than w. For example, we have
C[A] = 5 and C[B] = 8 by simple counting. The range
[C[w], C[w+1]) defines the suffix range R(w): R(A) = [5, 8),
for example (see that A appears as prefixes in [5, 8) in Fig. 2).
To understand how Algorithm 1 works, let us consider a
query P = BA. In Line 1, we have w = A, sp = 5, and
ep = 8. Consider the first (and last) iteration with i = 2. We
have sp = C[B] + 1 = 9 and ep = C[B] + 3 = 11 because
rankB (Tbwt , sp) = 1 and rankB (Tbwt , ep) = 3 by definition.
Therefore [sp, ep) = [9, 11) is returned at Line 7, which is
equivalent to R(BA) given in Fig. 2.
We can say that fast calculation of rankw enables the fast
execution of Algorithm 1 because all the operations except for

Input: BWT string of length n: Tbwt ,
Query string of length m: P
Output: Range of Tbwt that matches to P
w ← P [m − 1]; sp ← C[w]; ep ← C[w + 1]
for i ← 2 to m do
w ← P [m − i]
sp ← C[w] + rankw (Tbwt , sp)
ep ← C[w] + rankw (Tbwt , ep)
if sp ≥ ep then return NotFound
return [sp, ep)

rankw (Tbwt , i) are merely either substitutions or summations.
However, naı̈ve calculation of rankw with cumulative counting
incurs an unacceptable O(|Tbwt |) time.
4) Wavelet tree: A wavelet tree [10] storing Tbwt enables
fast calculation of rankw (Tbwt , i); its time complexity does not
depend on the data size |Tbwt |. Figure 3 illustrates a wavelet
tree for the string S = Tbwt in Eq. (2). The bit representation
of each symbol is predefined (e.g., Huffman coding based on
the frequency of each symbol in Tbwt ). Each node v in the
tree stores a bit vector Bv . For the root node v0 , Bv0 stores
the most significant bit (MSB) of each symbol in S. At the
second level, the symbols are divided into two parts based on
the bit value at the first level, while keeping the ordering. Each
bit vector stores the second MSB. Repeating such partitioning
recursively, we obtain the wavelet tree. In fact, Bv is stored
in a succinct dictionary [11], [12], which is a bit vector that
supports a bit-wise rank (i.e., rank0 (Bv , j) and rank1 (Bv , j))
in O(1) time.
There are several types of wavelet tree with different
compression characteristics that are determined by tree shape
and the type of succinct dictionary [13]. In CiNCT, we use a
Huffman-shaped wavelet tree (HWT) [14], whose tree shape
is that of the Huffman tree of S. It is known that an HWT can
compress a string S of length n to at most n(1+H0 (S))+o(n)
bits. Here, H0 (S) is the 0th order empirical entropy [15],
X nw
n
H0 (S) =
lg
,
(3)
n
nw
w∈Σ

where nw is the number of occurrences of w in S.
To calculate rankw (S, j), the wavelet tree calculates the
bit-wise rank value at each node v0 , v1 , · · · , vk between the
root and the leaf corresponding to the bit representation
w = b0 b1 · · · bk (see [10] for details). This indicates that bitwise rank operations required to obtain rankw (S, j) is equal
to k (i.e., the length of the bit representation of w). This fact
leads to the following result [13].
Theorem 1 (Rank on HWT): If rankw (S, j) is executed on
uniformly random w over S[0, n), it runs in O(1 + H0 (S))
time on average.
This result implies that a string with small entropy H0 (S)
achieves not only small size but also fast rank operation, which

$AAABDBBCCE$$$F#
1111000011011100

Level 1

Level 2

Level 3

Bit Repr.

{#,B,D,E,F}

{$,A,C}

BDBBEF#
0100111
{#,D,E,F}
{B}
DEF#
B
1001
{#,D}
{E,F}

$AAACC$$$
011111000

EF
01

Level 4
{E}

E

{A,C}

{$}

$

AAACC
11100
{A}
{C}

D#
10
{F} {#}

F

#

C

#:
$:
A:
B:
C:
D:
E:
F:

A

0110
10
111
00
110
0111
0100
0101

MSBs

{D}

D

Fig. 3. Wavelet tree: a bit representation of each symbol in a string S is
stored in a binary tree (this example is the HWT of the BWT of the trajectory
string Eq. (2)). Note that only bit vectors are stored in each node.

plays an important role in our theoretical analysis.
B. Compressed variants of FM-index
Let us consider a sub-path of length 3 in a real NCT dataset:
et−2 et−1 et . It is unlikely that two right turns occur in a row
because most vehicles go toward their destinations.
Considering such high-order correlations among symbols,
we can boost the compression. As noted before, the prefix
BA ∈Σ2 appears in [9, 11) (Fig. 2). The other prefixes W ∈
Σ2 have their corresponding ranges. Let us divide Tbwt based
on such prefixes W (called contexts of length two) as shown
in Fig. 4. These context blocks represent the next segment
et given the context W = et−1 et−2 . We have a chance of
compression because, as discussed above, the frequency of
symbols in each context is biased.
1) Compression boosting (CB): The above idea can be
generalized to any length of context. Let us divide Tbwt into
l blocks of context W ∈ Σk of length k: Tbwt = L1 L2 · · · Ll
(l ≤ σ k ). Storing each Lj in a 0-th order entropy compressor
such as an HWT, we can compress Tbwt to nHk (T ) + o(n).
Here, Hk is k-th order empirical entropy [15]:
X nW
Hk (T ) :=
H0 (TW ),
(4)
n
k
W ∈Σ

where TW is the concatenation of all symbols in T that precede
the context W . To support a fast rank operation on those
divided blocks, we need to precompute and store the rank
results at each location of l blocks for all w ∈ Σ.
Taking larger k seems to be desirable because Hk (T ) ≥
Hk+1 (T ) for all k ≥ 0 [15]. However, partitioning into many
blocks leads to the following problems in practice:
P1) Blocks of variable length lead to inefficient random
access to Tbwt .
P2) Index size increases because of the overhead of blockwise storage (e.g., pointers in Huffman trees).
P3) We have to save lσ integers for the rank results. This is
unrealistic for huge σ even if k = 1 (l = σ).
2) Variants of CB: There are some CB variants that
avoid the above problems. Fixed-block compression boosting
[16] adopts blocks of a fixed size. Although this solves P1
(and P2 partially), problem P3 remains for huge σ. Implicit
compression boosting (ICB) [17] avoids such explicit block
partitioning by using a compressed succinct dictionary called

Contexts
(length 2)

2nd column
1st column

last column (

)

F # CC D $$$ $ AA BB A B E
# $ $$ $ AAA B BB CC D E F
$ A AA B DBB C CE $$ $ F #

Fig. 4. Compression boosting of FM-index: Tbwt is divided into contexts
and each partition is compressed separately.

an RRR [12] in the wavelet tree of Tbwt . This implicit partition
solves P1 and P3. Brisaboa et al. [3] employed ICB to index
a trajectory string of NCTs. Specifically, they employed ICB
with a wavelet matrix [18], which is an efficient alternative to
a wavelet tree. We call this structure ICB-WM in this paper
(similarly, we refer to ICB with an HWT as ICB-Huff). As
discussed in our theoretical analysis, ICBs still suffer from
large overheads when applied to a string with large alphabet,
such as a trajectory string of NCTs.
III. P ROPOSED DATA STRUCTURE
A. Overview
For NCTs, the alphabet size σ can be millions because
it is the number of road segments in a road network. As
discussed in the previous section, this makes the compression
of trajectory strings inefficient, because the redundant bits in
wavelet trees increase as σ increases. To avoid this, we convert
trajectory strings into strings with a small alphabet via relative
movement labeling (RML), which is based on the sparsity of
road networks. Figure 5 and the following give an overview
of how to construct the proposed data structure, CiNCT.
1) Convert a set of NCTs into a trajectory string T .
2) Calculate the BWT of T and obtain Tbwt .
3) Construct an ET-graph GT and a relative movement
labeling (RML) function φ based on T (Section III-B)
4) Label Tbwt based on the RML function φ and obtain the
labeled BWT φ(Tbwt ) (Section III-C).
5) Store φ(Tbwt ) in an HWT with RRR and obtain the
proposed index structure (Section III-C).
As steps 1 and 2 are straightforward, we describe the details
of steps 3–5 in the following sections. We emphasize that the
NCTs are labeled after the BWT (step 4), otherwise we would
be unable to implement the suffix range query. Due to this
labeling step, we need to develop an algorithm that differs
from Algorithm 1. Such an algorithm is described in Section
IV. The theoretical consequences of CiNCT are described in
Section V. Here, we focus on the index structure.
Note that CiNCT basically deals with static data. We can
treat growing data by periodic reconstruction or by constructing an index for new data at certain time intervals.
B. Relative movement labeling (RML)
The RML converts trajectory strings into strings with small
alphabet based on the following fact: NCTs can only move between physically connected road segments. First, we describe
its idea based on the example in Fig. 1 (a). If a vehicle is on
a road segment w0 = A, the next segment w has to be B or
D. Hence, we label them 1 and 2, respectively. Generally, if
there are k connected road segments from a certain segment,

[INPUT]

NCTs

Step 1.

Existing methods

String
Step 2.

BWT

Labeled BWT

ET-graph

$
1

Step 3.

2
1

Labeling
function

B
A

Compressed Wavelet Tree

Fig. 5. Overview of CiNCT (index construction). Note that the existing ICBs
do not have the labeling step.

we can label them with 1, · · · , k. The sequences converted
with this relative movement labeling (RML) are expected to
have small alphabet because k is smaller than the maximum
out-degree of the road network. To define RML formally, let
us define an empirical transition graph (ET-graph).
Definition 3 (ET-graph): Let T be a string defined on an
alphabet Σ. An ET-graph GT of T is a directed graph that
satisfies: 1) the vertex set is Σ; 2) a directed edge (w0 , w) ∈
Σ × Σ exists iff there exists a substring ww0 in T . The edge
set is denoted by ET .
In other words, an edge exists iff a direct transition between
w0 and w exists in T . The ET-graph is a sparse graph
because it has a similar topology to the original road network.
Figure 6 (a) illustrates the ET-graph of the trajectory string
T given in Eq. (1). Note that ET-graphs include the special
symbols $ and #.
1) Definition of RML: The RML can be defined as an
integer assigned on each edge of ET-graph (see Fig. 6 (a)).
For example, the transition A → B is labeled 1. The transition A → D must have the different label, otherwise we
cannot distinguish them. For transition w0 → w, we denote
such a labeling function by φ(w|w0 ). For example, we have
φ(B|A) = 1 and φ(D|A) = 2. To make the labeling distinct
based on the previous symbol w0 , the RML function φ must
satisfy the following requirement.
0
• Requirement: The RML function φ(·|w ) must be a one0
to-one map for any w .
Now, we discuss how to construct the RML function φ that
satisfies the requirement above. Let us consider the out-vertex
set of w0 , defined as Nout (w0 ) = {w|(w0 , w) ∈ ET }, that
determines the set of vertexes directly accessible from w0 .
Based on the ET-graph and out-vertex set, we define φ(·|w0 )
as follows. Given w0 , assign a different small integer cww0
to each w ∈ Nout (w0 ) and define φ(w|w0 ) := cww0 . It is
clear that φ(·|w0 ) is a one-to-one map. If w ∈
/ Nout (w0 ), we
0
cannot define φ(w|w ). However, this is not a problem because
w ∈
/ Nout (w0 ) indicates that the string ww0 is not found in
T , which tells us the result of pattern matching is null. This
point is important for our search algorithm.
2) Finding an optimal RML: The RML φ described above
does not define a unique labeling function because we have
not yet specified a concrete way to assign the small integers
cww0 . Here, we propose a strategy based on a bigram count
nww0 (i.e., the frequency of ww0 in T ). The elements in
Nout (w0 ) are sorted in descending order of bigrams nww0 . The
vertex w with the largest bigram count is given the smallest

#

1

1
Context
(k=1)

C

# $$$$ AAA BBB CC D E F

2

1

Step 5.

1
1

D

2

F

1

E

(a) ET-graph and RML

$ AAAB DBB CCE $$ $ F #
1 1112 211 112 11 1 1 1

(b) Labeling Tbwt with RML

Fig. 6. (a) ET-graph of our example Eq. (1): each node represents a road
segment, and an edge exists if the corresponding transition occurs in T .
The integer on each edge is the corresponding label. (b) Relative movement
labeling: φ(Tbwt ) produces a string with a lower entropy than that of Tbwt .

label, 1. The second-most frequent vertex is labeled 2, the
third-most frequent vertex is labeled 3, and so on. The labels
shown in Fig. 6 (a) are determined in this way. For example,
since we have nBA > nDA (nBA = 2 and nDA = 1), the
edge from A to B has the smallest label 1: φ(B|A) = 1.
Applying a labeling scheme shown in the next section, this
labeling strategy generates a low-entropy sequence φ(Tbwt ) as
shown in Fig. 6 (b), because the distribution of the resulting
symbols is biased toward smaller integers (i.e., 1 is the largest
fraction). For this example, we have H0 (Tbwt ) = 2.8 and
H0 (φ(Tbwt )) = 0.7 (unit: bits).
One might wonder whether there exists a better labeling
strategy. We prove, however, the optimality of the labeling that
leads to strong conclusions: our RML achieves the smallest
size and the fastest search. See Section V-A for details.
C. Data structure
Here, we describe how to obtain φ(Tbwt ) and the final index
(steps 4 and 5 in Section III-A).
1) Labeling BWT (step 4): Based on the RML function φ
obtained in the previous section, the BWT Tbwt is converted
to φ(Tbwt ) in the following manner. For example, let us focus
on the third block of Tbwt , DBB, in Fig. 6 (b). This block
corresponds to the context of A, which indicates that the
previous symbol of these DBB is A. Hence, DBB is labeled
as 211 because φ(B|A) = 1 and φ(D|A) = 2 in Fig. 6 (a).
All the other blocks also can be labeled in the same manner.
2) Storing to a compressed wavelet tree (step 5): In this
step, we store the labeled BWT φ(Tbwt ) to an HWT. For
bit vectors in an HWT, we adopt a practical version of
the compressed succinct dictionary called RRR [19]. This is
a straightforward step. Figure 7 depicts the comparison of
Huffman trees of Tbwt and φ(Tbwt ) for the example in Figure
6 (b). The Huffman tree of φ(Tbwt ) is obviously simpler than
that of Tbwt . Because these tree shapes are the same as those
of HWTs, this simplification explains intuitively why CiNCT
is small and fast. For more details, see Section V.
An RRR bit vector has one parameter b, that controls the
size of the internal blocks. For larger b, we obtain better
compression but slower search (rank calculation) in general,
and vice versa. This b is the only parameter in CiNCT.
However, in Section VI, we show that this parameter has only
a small influence on the index size and the search time.

Bit Repr.
#:
$:
A:
B:
C:
D:
E:
F:

0110
10
111
00
110
0111
0100
0101

0
0

1

1
0

B

1

0
1 $

0

1 0

1

E

F #

D

(a) Tbwt

Bit Repr.

0

1

C

A

1: 0
2: 1

0

$ AAAB DBB CCE $$ $ F #

1

1

One to one correspondence
in
(shaded area)

2

1 1112 211 112 11 1 1 1
(b) φ(Tbwt ): CiNCT

Fig. 7.
Huffman trees of Tbwt and φ(Tbwt ): each leaf corresponds to
a symbol. The tree generated by CiNCT (b) is simpler than that from the
existing technique (a).

3) Storing ET-graph: We use an adjacency list to represent
the ET-graph GT . The value φ(w|w0 ) is assigned to the edge
(w0 , w) ∈ ET . Thus we can obtain φ(w|w0 ) in O(δ) time by
a linear search over Nout (w0 ). We also assign C[w] to each
vertex w in GT . Correction terms Zw0 w , which are introduced
in Section IV-A, are also attached to (w0 , w) ∈ ET . Note that,
since GT is sparse, the space needed to store GT is negligible
when |T | gets large.
IV. P ROPOSED QUERY PROCESSING ALGORITHMS
Here, we describe another key concept of this paper, PseudoRank, then show algorithms for two types of queries, suffix
range queries and sub-path extraction queries.
A. PseudoRank
As mentioned in Section II-A3, fast calculation of
rankw (Tbwt , j) is needed for Algorithm 1. The original FMindex stores Tbwt in a wavelet tree to calculate ranks quickly.
In our case, however, we do not have the original Tbwt but only
have the labeled φ(Tbwt ). Can we obtain the rank values for
the original BWT by using only the labeled BWT? Seemingly,
this is difficult because different symbols are mapped to the
same label (e.g., both A and C are converted to 1 as illustrated
in Fig. 6 (b)).
The key idea in CiNCT is to simulate the rank operation
over Tbwt . Figure 8 illustrates this idea. Let us consider the
range R(A) = [C[A], C[B]) and j ∈ R(A). Because the
substring Tbwt [C[A], C[B]) = DBB is labeled as 211 by using
the one-to-one map φ(·|A) as described in Section III-C, the
following two counts are equivalent for ∀j ∈ R(A):
0
• the number of occurrences of D within the range R :=
[C[A], j) in Tbwt (the shaded region in Fig. 8), and
0
• the number of occurrences of 2 within R in φ(Tbwt ).
This balancing relationship holds in general. Let us consider
a context w0 . For all j such that C[w0 ] ≤ j ≤ C[w0 + 1],
let us consider a range R0 := [C[w0 ], j). For a symbol w ∈
Nout (w0 ), the number of occurrences w within R0 in Tbwt and
that of the label η := φ(w|w0 ) within R0 in φ(Tbwt ) are the
same because of the one-to-one requirement for φ(·|w0 ). This
leads to the following balancing equation:
rankw (Tbwt , j) − rankw (Tbwt , C[w0 ])
= rankη (φ(Tbwt ), j) − rankη (φ(Tbwt ), C[w0 ]).

(5)

Fig. 8.

Basis of the balancing equation (Eq. (5)) for PseudoRank

Rearranging this equation, we have the following theorem,
which allows us to simulate the rank operation.
Theorem 2 (Pseudo-rank): If w ∈ Nout (w0 ) and
C[w0 ] ≤ j ≤ C[w0 + 1], then we have
rankw (Tbwt , j) = rankη (φ(Tbwt ), j) − Zw0 w ,

(6)

0

where η := φ(w|w ) and
Zw0 w := rankη (φ(Tbwt ), C[w0 ]) − rankw (Tbwt , C[w0 ]). (7)
We emphasize that the correction term Zw0 w does not
depend on j, implying that the number of correction terms
needed is equal to |ET |. Importantly, this property allows us to
precompute and store the correction terms (as noted in Section
III-C, they are attached to each edge (w0 , w) ∈ ET ).
This theorem produces Algorithm 2, which calculates the
rank values using only φ(Tbwt ). We also emphasize that
PseudoRank does not allow us to calculate rank values for
all pairs of (w, j). However, this limitation is not a problem
for our search algorithm, as shown in the next subsection.
Algorithm 2: Pseudo calculation of rankw (Tbwt , j) by using
only φ(Tbwt ) (PseudoRank(φ(Tbwt ), j, w, w0 , Zw0 w ))

1
2
3
4

Input: Labeled BWT string of length n: φ(Tbwt ),
Location of rank j, Correction term Zw0 w ,
Target symbol w, Previous symbol w0
Output: The value of rankw (Tbwt , j)
if w ∈ Nout (w0 ) and C[w0 ] ≤ j ≤ C[w0 + 1] then
η ← φ(w|w0 )
// RML
return rankη (φ(Tbwt ), j) − Zw0 w
return NotFound

B. Suffix range query with CiNCT
With the PseudoRank, we can simulate rankw (Tbwt , j) using
only the wavelet tree of φ(Tbwt ) and the correction term Zw0 w
(Eq. (7)). Replacing the rank operations in Algorithm 1 with
PseudoRank, we obtain our search algorithm (Algorithm 3),
whose correctness is shown below.
1) Correctness of the algorithm: To guarantee that Algorithm 3 is equivalent to Algorithm 1, we have to check the
following two conditions on PseudoRank (Theorem 2) are
satisfied immediately before Line 7: (c1) w ∈ Nout (w0 ); (c2)
C[w0 ] ≤ sp ≤ C[w0 + 1] and C[w0 ] ≤ ep ≤ C[w0 + 1].
As noted previously, no substring ww0 appears in T if w ∈
/
Nout (w0 ); hence, NotFound is returned if w ∈
/ Nout (w0 )

Algorithm 3: Finding the suffix range [sp, ep) for a given query

Algorithm 4: Extracting a sub-path T [i − l + 1] · · · T [i − 1]T [i]
for given j = ISA[i] and l > 0 (extract)

P of length m based on φ(Tbwt ) (LabeledSearchFM)

1
2
3
4
5
6
7
8
9
10
11

Input: Labeled BWT string of length n: φ(Tbwt ),
Query of length m: P [0, m),
Correction terms: {Zw0 w }
Output: Range of Tbwt that matches to P
w ← P [m − 1]; sp ← C[w]; ep ← C[w + 1]
for i ← 2 to m do
w0 ← w
// Save previous symbol
w ← P [m − i]
if w ∈
/ Nout (w0 ) then
return NotFound
sp ← C[w]+PseudoRank(φ(Tbwt ), sp, w, w0 , Zw0 w )
ep ← C[w]+PseudoRank(φ(Tbwt ), ep, w, w0 , Zw0 w )
if sp ≥ ep then
return NotFound
return [sp, ep)

at Line 6. Therefore, (c1) w ∈ Nout (w0 ) holds immediately
before Line 7. For (c2), before Line 7, sp satisfies
sp = C[w0 ] + rankw0 (Tbwt , sp0 ),

(8)

where sp0 is the previous value. By the rank definition,
0 ≤ rankw0 (Tbwt , j) ≤ C[w0 + 1] − C[w0 ],

0 ≤ ∀j < |T |

holds, where C[w0 + 1] − C[w0 ] means the number of occurrences of w0 in T . Combining this inequality with Eq. (8), we
obtain C[w0 ] ≤ sp ≤ C[w0 + 1]. We can prove the condition
for ep in a similar manner.

1
2
3
4
5
6

Input: Labeled BWT: φ(Tbwt ), Position on Tbwt : j,
Extraction length: l, Correction terms: {Zw0 w }
Output: A substring S := T [i − l + 1] · · · T [i − 1]T [i]
w0 ← BinarySearch(j, {C[w0 ]})
// T[i]
for k ← 1 to l do
η ← φ(Tbwt )[j]; w ← decode(η|w0 ); S[l − k] ← w
j ← C[w] + PseudoRank(φ(Tbwt ), j, w, w0 , Zw0 w )
w0 ← w
// Save previous symbol
return S

V. T HEORETICAL ANALYSIS
In this section, we explain theoretically why CiNCT is
compact and fast. We first show the optimality of our proposed
RML, that is, the labeled BWT φ(Tbwt ) achieves the smallest
entropy. Then, we explain that such a small entropy contributes
high compressibility and fast query processing. We also show
that RML is better than other labeling method called MEL,
recently proposed in [1].
A. Optimality of RML
The 0th order empirical entropy H0 given in Eq. (3)
plays important roles in our analysis. First, we show the
labeling strategy based on bigram counts nww0 proposed in
Section III-B achieves the minimum value of H0 among all
possible labelings.
Theorem 3 (Optimality): Let φ∗ be the RML based on the
bigram ordering strategy and φ be any possible RML that
satisfies the requirement in Section III-B. Then, we have
H0 (φ∗ (Tbwt )) ≤ H0 (φ(Tbwt )).

(9)

C. Extracting a sub-path with CiNCT
Here, we describe another important query, sub-path extraction query. For example, let us focus on the third sorted
rotation in Fig. 2. Its suffix of length four is FEBA (colored in
blue). This corresponds to the example NCT T1r in Eq. (1). In
this way, the sub-path extraction queries recover a sub-path of
length l from an arbitrary position j in BWT Tbwt (j = 3 for
the example above). This query is useful if we need to obtain
certain NCTs stored in BWT string, or we need to recover the
entire trajectory string. Formally, extract(j, l) returns T [i−l, i)
where i = SA[j] (SA is the suffix array of T ). The subscript
j is often referred to as inverse suffix array (j = ISA[i]).
Algorithm 4 shows how to obtain extract(j, l) using only
φ(Tbwt ) and the ET-graph. This is obtained by mimicking
LF-mapping [8] with PseudoRank. Line 1 performs a binary search to find the last character T [i] = w0 such that
C[w0 ] ≤ j < C[w0 + 1]. Line 4 first accesses the j-th
character of φ(Tbwt ) (i.e., the labeled Tbwt [j]), then decodes
the Tbwt [j] = T [i − k − 1] = w using the ET-graph. Line 5
is similar to Line 7 in Algorithm 3, which jumps to the next
position on Tbwt (LF-mapping simulated by PseudoRank).

Proof: See Appendix A.
As a special case of this theorem, we obtain an unlabeled
case result, i.e., H0 (φ∗ (Tbwt )) ≤ H0 (Tbwt ), by putting as
φ = id (identity labeling). Importantly, we see that
H0 (φ∗ (Tbwt ))  H0 (Tbwt )

(10)

holds for real NCT datasets in our experiments (Table III).
B. Compressed size
1) Evaluating space overheads: The data structure of
CiNCT consists of two parts: the labeled BWT φ(Tbwt ) and
the ET-graph GT . As noted in Section III-C, the size of GT
is negligible when |T | is large. Here, we compare the sizes of
Tbwt and φ(Tbwt ) stored in HWTs with RRR. Note that these
corresponds ICB-Huff and CiNCT, respectively. The main
advantage of CiNCT comes from the lower space overhead
due to RRR, as explained below. For a given bit vector B, the
practical RRR with the parameter b uses at most
|B|H0 (B) + |B| · h(b)

(11)

[19]. We call the second term the
bits 1 where h(b) = lg(b+1)
b
RRR-overhead. For b = 63, we have an overhead of h(b) =
(lg 64)/63 ' 0.095 bits per bit.
For a given string S, the average code length with Huffman
coding is at most (1 + H0 (S)) P
bits. Hence, the total length
of bit vectors in the HWT is
v |Bv | ' |S|(1 + H0 (S))
(Section II-A). Summing the RRR-overheads over all internal
nodes v in the HWT, we obtain total bits of the overhead:
X
|Bv | · h(b) ' |S|(1 + H0 (S)) · h(b).
(12)
v

The right-hand side implies that the RRR-overhead of a
sequence S is small if its entropy H0 (S) is small. Therefore,
Eq. (10) indicates that the space overhead for CiNCT is much
smaller than that for ICB-Huff.
2) High-order compression: Here, we analyze the remaining first (and dominant) term in Eq. (11). Summing this term
over all internal nodes v in the HWT, we find that the total bits
needed for this term achieves k-th order entropy Eq. (4) for
all k > 0. This property implies that our method guarantees a
high compressibility in information theoretic sense. Note that
this kind of entropic bound has not been guaranteed by the
existing shortest-path based NCT compressors.
Theorem 4: For all k > 0, the total bits required to store
φ(Tbwt ) in an HWT with RRR, apart from the overhead
Eq.(12), are |T |Hk (T ) + O(lσb), where l ≤ σ k is the number
of distinct contexts W ∈ Σk in T .
Proof: See Appendix B.
C. Processing time of suffix range queries
To evaluate whether Algorithm 3 is faster than Algorithm 1,
we focus on the time complexity of the rank operation. As
stated in Theorem 1, rankw (S, j) runs in O(1 + H0 (S))
time2 . Hence, the relationship H0 (φ(Tbwt ))  H0 (Tbwt )
again explains why CiNCT is faster than ICB-Huff. Of course,
Algorithm 3 incurs an additional cost in calculating φ(w|w0 ),
but this is not serious for a sparse GT .
Moreover, we have the following theorem implying that
the search time does not depend on the road network size
σ but depends only on the maximum out-degree δ of the road
network (which is usually less than four).
Theorem 5 (σ-independence): Let P ∈ E ∗ be any query
path ($ is not included). Algorithm 3 runs in O(|P | · δb) time.
Proof: For any w, w0 ∈ E, we have η := φ(w|w0 ) ≤ δ+2.
By the construction of RML, η is at least the δ + 2-th most
frequent symbol in φ(Tbwt ). Thus η is at most located at the
δ + 2 level of the Huffman tree. Hence, rankη (φ(Tbwt ), j)
in Eq. (6) runs in O(δb) time (remember the bit-wise rank
operation in practical RRR [19] requires O(b) time). Since
PseudoRank is calculated at most 2|P | − 2 times in Algorithm
3, this leads to the conclusion.
Other FM-indexes do not satisfy this property. Note that this
time complexity also does not depend on the data size |T |.

D. Comparison of RML with MEL
Minimum entropy labeling (MEL) is a labeling scheme for
NCTs that was recently proposed in [1], which works as a
preprocessor for general compressors, such as Huffman coding
or LZ coding (i.e., pattern matching was not considered).
Similar to RML, MEL converts a sequence of road edges to
a low entropy sequence of small integers as follows:
w1 w2 · · · wn → ψ(w1 )ψ(w2 ) · · · ψ(wn )

where ψ : E → N is the MEL function. Different labels are
assigned to road segments that shares head node v (Fig. 9(b)).
By contrast, our RML conversion is as follows:
w1 w2 · · · wn → φ(w1 |$)φ(w2 |w1 ) · · · φ(wn |wn−1 ).

(14)

Unlike RML, the MEL function ψ does not consider the previous symbol. Specifically, MEL labels based on the unigram
frequencies, which are shown as nA and nB in Fig. 9(b).
Conversely, our RML, shown in Fig. 9(a), is based on bigram
frequencies, nXA , nXB , nYA , and nYB .
Given these differences, the advantage of RML can be
intuitively explained as follows. Real trajectories tend to go
straight rather than turn left or right, as shown in Fig. 9.
Because RML considers the previous road segment, RML can
take account the direction of the movement, whereas such
information is lost in MEL. This implies that RML can capture
a higher-order correlation compared to MEL. Although MEL
also has the optimality of entropy, it cannot be better than
RML. The experimental comparison is shown in Section VI-D.
Mathematically, we have the following theorem.
Theorem 6: For any trajectory string T , RML achieves a
smaller 0th order empirical entropy than MEL does.
Proof: Considering the size of the feasible labeling space,
we find that our labeling space {φ(w|w0 )} is a superset of
that of MEL, {ψ(w)}. In other words, MEL can be emulated
by an RML φ̄ that is not necessarily optimal. Therefore, the
optimality of RML (Theorem 3) leads to the conclusion.
VI. E XPERIMENTS
A. Experimental setup
1) Implementation: All methods were implemented in C++
and compiled with g++ (version 4.8.4) with the -O3 option. We
used the sdsl-lite library (version 2.0.1) for (in-memory)
wavelet trees (http://github.com/simongog/sdsl-lite/). The BWT
was calculated using sais.hxx (http://sites.google.com/site/
yuta256/sais/). Experiments were conducted on a workstation
with the following specifications: Intel Core i7-K5930 3.5GHz
CPU (64-bit, 12 cores, L1 64kB×12, L2 256kB×12, L3
15MB), DDR4 32GB RAM, Ubuntu Linux 14.04.
(a)
A

Y

nXA

nYA

X

B

nYB

1 In

fact, there exist non-dominant terms that are not included in
this equation. See [19] for details.
2 To be exact, this complexity is proportional to b because practical
RRR [19] runs bit-wise rank in O(b) time.

(13)

nXB

(b)

nA=nXA+nYA

A
v
B

nB=nXB+nYB

Fig. 9.

Comparing two NCT labeling methods: (a) RML; and (b) MEL [1]

2) Competitors: Table II lists the competitors used in this
paper. We used five FM-index variants: uncompressed (UFMI,
FM-GMR) and compressed (ICB-WM, ICB-Huff, FM-APHYB). The block-size parameter b had to be specified for
CiNCT, ICB-Huff, and ICB-WM. Unless otherwise noted,
we use b = 63. FM-GMR [20] and FM-AP-HYB [21]
are FM-index variants that are tailored for huge σ and that
support O(log log σ) rank operation (faster than the O(log σ)
of UFMI); they are available in sdsl-lite library. These
were the fastest (FM-GMR) and the smallest (FM-AP-HYB)
methods for huge σ in a recent benchmark [22].
There are many possibilities for compressing NCTs by
combining simple techniques such as run-length encoding.
However, we do not consider such techniques in this study
because pattern matching is not supported in sublinear time.
In our prior evaluation, the Boyer-Moore method (linear time
search) was at least four orders of magnitude slower than
CiNCT even if T was stored in an in-memory uncompressed
array. In this study, we thus only consider RePair [23], a
standard benchmark in stringology which showed the best
compression ratio in the initial evaluation, and PRESS [24],
which is the shortest-path-based NCT compressor (Note that
‘Q?’ is unchecked in Table II for these methods).
3) Measurement: The search time was averaged over 500
suffix range queries of length 20 randomly sampled from the
data. For evaluation of the data size of CiNCT, the size of the
ET-graph is included.
4) Datasets: The datasets used in this study are as follows:
• Singapore: NCTs of taxi cabs used in [24]. This dataset
contains many transitions without physical connection.
• Singapore-2: Preprocessed Singapore dataset such that
transitions between two road segments without a physical
connection are interpolated with the shortest path. This
is used to evaluate the gain against the noisy dataset
(Singapore).
• Roma: GPS trajectories of taxi cabs in Rome. NCT
representations were obtained by HMM map-matching
[25] (http://crawdad.org/roma/taxi/).
• MO-gen: NCTs generated by the moving object generator (http://iapg.jade-hs.de/personen/brinkhoff/generator/).
• Chess: All chess game records (Blitz, 2006–2015, 1.87
million games, http://www.ficsgames.org). Each openings
(10 moves) is converted into a string of hash values of
Forsyth-Edwards notation.
Although Chess is not a vehicular dataset, it also has a
sparse ET-graph GT because of the characteristics of chess
games. This is included to show the possibility that CiNCT
is applicable to targets other than NCTs. Table III lists the
statistics of the datasets, which are used to explain the results.
B. Comparison with various FM-indexes
Evaluation results for data size and processing time of
suffix range queries are shown in Fig. 10. We observe that
CiNCT requires less than 2 bits per symbol to store NCTs,
and pattern matching of length 20 is processed in a few tens of

TABLE II
O UR PROPOSED METHOD AND ITS COMPETITORS∗
Method
CiNCT
UFMI

Data
φ(Tbwt )
Tbwt

ICB-WM
ICB-Huff
FM-GMR

Tbwt
Tbwt
Tbwt

FM-APHYB
PRESS
[24]
Re-Pair
[23]

Tbwt

∗
†


T
T

Description
HWT with RRR
WM [18] with uncompressed bitmap [11]
WM with RRR [18]
HWT with RRR [17]
FM-index for huge σ with
O(log log σ) rank [20]
FM-index for huge σ with
O(log log σ) rank [21]
The state-of-the-art trajectory compressor
A standard benchmark compressor in stringology

C?†

Q?‡

For the first four method, the type of WT used is in description
Uncompressed or compressed / ‡ Supports suffix range query or not
WM: wavelet matrix

TABLE III
S TATISTICS OF EACH DATASET
Dataset
Singapore
Singapore-2
Roma
MO-Gen
Chess
†
‡

|T |
53M
75M
12M
193M
20M

lg σ
15.5
15.5
15.5
17.4
18.8

H0 (T )
13.8
14.0
13.0
13.0
10.3

H0 (φ)†
1.8
1.3
0.9
2.8
2.0

H1 (T )
1.5
1.1
0.7
2.5
1.4

d¯ ‡
26.8
4.0
2.4
8.8
1.6

H0 (φ) means H0 (φ(Tbwt ))
d¯ is the average out-degree of the ET-graph GT .

microseconds. We also observe that CiNCT outperforms the
competitors in terms of both data size and query processing
time. We explain these results in detail below.
1) Data size: Compared with ICB-Huff and ICB-WM,
CiNCT reduces the data size by up to 78% and 57%, respectively. As explained in Section V, the space overhead
decreases if H0 (S) decreases. From Table III we can confirm
that H0 (φ(Tbwt ))  H0 (Tbwt ) holds for all datasets (note
that H0 (T ) = H0 (Tbwt )). This explains why CiNCT shows
this significant improvement.CiNCT even shows better compression than the smallest variant FM-AP-HYB, which was
designed for huge σ. The improvement in Singapore-2 is
larger than that of Singapore. Because “gapped” transitions
are interpolated in Singapore-2, the ET-graph gets sparser
(d¯ = 26.8 → 4 in Table III). This reduces the overhead regarding the ET-graph (this is confirmed through the difference of
CiNCT and CiNCT (w/o ET-graph)).
2) Processing time of suffix range queries: CiNCT is
always much faster than ICB-Huff and ICB-WM; the speedups
are up to 7 and 25 times, respectively. Surprisingly, CiNCT is
even faster than those of the uncompressed indexes (UFMI
and FM-GMR). Again, this speedup can be explained by
the shallowness of the HWT of CiNCT. Because we have
H0 (φ(Tbwt ))  H0 (Tbwt ), the HWT becomes shallower. This
decreases the number of bit-wise rank operations in the HWT,
as pointed out in Section V-C.
3) Effect of block size b: As mentioned in Section III-C,
as b becomes larger, the results show better compression but
slower search. However, as shown in Fig. 10, the sensitivity
to the block size parameter b is very small for CiNCT. This

ICB-Huﬀ
b=63
b=31

b=63

64

b=31

b=15

FM-AP-HYB
b=63
b=31
b=15

b=15

b=15

64

b=15

1

2

4

8

16

32

64

1

Space (bits/symbol)

2

4

8

16

32

64

(c) Roma
0.5

Space (bits/symbol)

1

2

4

8

16

32

(d) MO-gen
1

Space (bits/symbol)

2

4

8

16

32

(e) Chess

1

4

(b) Singapore-2

16

4

16

(a) Singapore

4

16

16

64

64

FM-GMR
b=63
b=31

16

256

256

b=31

4

Search time ( µs)

b=63

ICB-WM

64

UFMI

256

CiNCT (w/o ET-graph)

256

CiNCT

64

1

2

Space (bits/symbol)

4

8

16

32

64

Space (bits/symbol)

0

5

10

15

20

Fig. 11. |P | vs. search time:
(Singapore dataset)

15

2 2 16 2 17 2 18
Alphabet size: σ

Fig. 12. CiNCT shows the best σ-dependence (Left:
index size, right: search time / RandWalk dataset)

indicates that the proposed method is nearly parameter-free.
4) Effect of |P |: Figure 11 shows how the processing
time of suffix range queries increases as the query length
|P | increases. For all methods, the processing time grows
linearly as the query length increases, because the numbers of
iterations in Algorithm 1 and Algorithm 3 is proportional to
|P |. We observe that CiNCT shows the slowest growth among
all methods.

Table IV compares the compression ratio, which is defined
as the uncompressed size (binary file of 32-bit integers)
divided by the compressed size. We observe that CiNCT shows
better compression than the existing methods. In particular,
our method is better than MEL, which also showed the best
compressibility in recent evaluation of NCT compression [1].
Note that the road network storage is not included in MEL
evaluations whereas it is considered for CiNCT (as ET-graph).
D. Effect of labeling strategy
1) Comparison with MEL: According to our analysis in
Section V-D, RML achieves lower entropy than MEL does.
We show our experimental results from two real NCT datasets,
i.e., Singapore2 and Roma. Table V provides a comparison
of the entropy achieved by RML and MEL. These results show
that our RML obtained approximately a 30% smaller entropy
than that of MEL.
2) Optimality: In Section III-B, we proposed a labeling
strategy that assigns small integers cww0 sorted by the bigram
counts nww0 . The data size and search time under this strategy

28
27
26

UFMI
FM−GMR
CiNCT
FM−AP−HYB

23

CiNCT (w/o ET-graph)

25

Search time ( µs)

CiNCT
FM−AP−HYB

ICB−Huﬀ

24

22
21
20

ICB−Huﬀ
ICB−WM

ICB−WM

22 23 24 25 26
Avg out−degree d

22 23 24 25 26
Avg out−degree d

Fig. 13. Dependence on out-degree (Left: index size,
right: search time / RandWalk dataset)

TABLE IV
C OMPRESSION RATIO ( LARGER IS BETTER )
CiNCT
MEL†
Re-Pair
bzip2
PRESS‡
zip
†
‡

C. Comparison with several compression methods

FM−GMR

2 −1

2
2

5

2 14 2 15 2 16 2 17 2 18
Alphabet size: σ

UFMI

2 −2

2

7
6

UFMI
CiNCT
FM−GMR

CiNCT (w/o)
14

FM−AP−HYB

Index size (bits / symbol)

FM−AP−HYB

ICB−WM

2

CiNCT

ICB−Huﬀ

4

ICB−WM

ICB−Huﬀ

Search time ( µs)

FM−GMR

2

Query length: |P|

2

8

Index size (bits / symbol)
2
2 −1
20
21
22

UFMI

−2

50

100

150

CiNCT
UFMI
ICB−Huﬀ
ICB−WM
FM−GMR
FM−AP−HYB

0

Search time ( µs)

200

Fig. 10. Data size/search time (suffix range query): the proposed method shows the best performance. CiNCT (w/o ET-graph) is used to show the data size
without ET-graph. The block size used in the RRR bit vectors is parameterized as b ∈ {15, 31, 63}. Results for the other methods in Table II were omitted
because their linear-time search was too slow.

Singapore
10.5
N/A
8.4
5.3
4.6
2.5

Singapore-2
27.0
15.8
11.4
5.6
N/A
2.5

Roma
25.2
21.2
20.6
13.6
N/A
5.0

Mo-Gen
25.6
N/A
20.6
5.3
N/A
2.6

Chess
10.3
N/A
11.0
7.1
N/A
3.9

Huffman coding was used after labeling, as in [1]. We evaluated only for ungapped
datasets because MEL assumes no gap (see Singapore-2 explanation in Sec.??).
Only the result for the Singapore dataset [24] is shown because no available
implementation was found.

TABLE V
C OMPARISON OF ENTROPY (RML AND MEL)
Dataset
Singapore2
Roma

RML (Proposed)
1.26
0.76

MEL [1]
1.93
0.99

are expected to be better than those of any other possible
labeling strategy, because we showed the optimality of our
strategy (Theorem 3). Here, we compare our strategy with
the random sorting strategy, which assigns randomly shuffled
small integers cww0 ∈ {1, · · · , |Nout (w0 )|}. Figure 14 shows
the comparison for the five datasets (b ∈ {15, 31, 63}). We
observe that the index size and the search time of the bigram
sorting strategy are always better than those of random sorting
strategy. Compared to the random strategy, it reduces the data
size by up to 32%, and the search time by up to 57%. These
results indicate the importance of the bigram sorting strategy.

Fig. 15.

Extraction time

E. Effect of ET-graph size/shape
1) Effect of σ: In Theorem 5, we showed that the search
time of CiNCT does not depend on the size σ of the road map.
Here, we investigate what happens when σ grows. For the
experiment, we use synthetic data RandWalk: random walks
on a directed random Poisson graph. The average out-degree
d¯ of the graphs is fixed at four, and |T | is set to 800σ.
In Fig. 12, CiNCT shows good scalability against σ,
whereas the index sizes and the search times of the existing
methods both increase. The search time of CiNCT is almost
constant, as predicted by Theorem 5.3 The other methods do
not show this property. For example, both the index size and
the search time of UFMI at σ = 218 are 30% larger compared
to the σ = 214 case.
2) Effect of sparsity: Next, we investigate how the data
¯
size and search time behave against the average out-degree d.
Figure 13 shows the results for the RandWalk dataset used in
Section VI-E. For comparison, we fixed σ = 216 and |T | =
100M, and changed d¯ between 22 and 27 . We observe that the
sparsity of the ET-graph is the key factor for CiNCT. Although
the compression performance of CiNCT is the best, the data
size grows quickly. This is due to two factors: the increase
of ET-graph size and the increase of the depth of HWT. This
result is a natural consequence of our assumption that the road
network is highly sparse. However, this result shows that our
method works for larger d¯ than in the road network case, d¯ '
22 . This result opens the door to applications to datasets not
mentioned in this paper (e.g., symbol-valued time series).
F. Sub-path extraction time
Here, we evaluate extract queries described in Section IV-C.
We evaluated the extraction time for obtaining the entire T ,
that is, l = |T | and j = 0. Figure 15 compares the extraction
times (per symbol) for the four datasets. We observe that
CiNCT shows the fastest extraction among the competitors
(twice as fast as UFMI). Again, this can be explained by the
fast rank calculation in CiNCT (PseudoRank), as discussed
above. Note that we omitted the results for FM-AP-HYB
because random access to Tbwt was not supported in the
sdsl-lite library.
3 In fact, the search time of CiNCT increases slightly because |T |
increases with σ in our setting, leading to a lower cache hit ratio.
We confirmed the exact constant search time when |T | is constant.

80
60
40
20

FM−AP−HYB

UFMI

CiNCT

FM−GMR

Comparison of labeling strategies

ICB−WM

Fig. 14.

ICB−Huff

23

ET−graph−build
WT−build
BWT

0

22

Space (bits/symbol)

Index construction time (s)

21

Chess

Chess

20

b=31
b=15

ICB−Huff
ICB−WM

MO−gen

Bigram sorting (Proposed)
Random sorting

CiNCT
UFMI
FM−GMR

Roma

b=63

Roma

b=15

Singapore

b=31

12

MO−gen

Singapore

8

b=15

4

b=31

0

b=63

b=63

Extraction time (µs/symbol)

Search time ( µs)

20 21 22 23 24 25 26

Singapore-2

Fig. 16. Index construction time
(Singapore dataset)

G. Index construction time
Figure 16 compares the index construction times of FMindex variants. The construction time of CiNCT is comparable
to that of ICB-Huff, and shorter than those of the other
methods. In fact, CiNCT is the second fastest among the
considered competitors. ET-graph-build in Fig. 16 includes
all operations that are not needed for the other methods. Here,
we can see the overhead for the construction of the ET-graph
is not a serious problem. Note that all additional operations,
including the construction of GT from T , obtaining RML
function φ, labeling Tbwt , and calculation of Zw0 w , can be
executed in linear time O(|T |), which implies the scalability
of construction.
VII. R ELATED WORK
1) Trajectory indexing: One important application of our
method is trajectory indexing. Although there are numerous
studies on this topic as shown in a survey [26], we present only
the most relevant ones here. MON-tree [27] is one of the most
famous methods to index NCTs. This method, as well as many
of other NCT indexing methods, mainly focuses on spatiotemporal range queries. Krogh et al. [28] proposed a data
structure similar to MON-tree designed for a different type of
query called strict path query, which aims to find trajectories
in a database that traveled along a given sub-path P during a
given time interval I. This method was implemented based on
B+ -trees. Koide et al. [6] showed that strict path queries can be
efficiently processed using suffix range queries. This method
is a hybrid data structure that indexes timestamps and spatial
paths using B+ -trees and FM-index, respectively. Brisaboa et
al. [3] also employed a compressed suffix array to store spatial
paths. As for timestamps, an in-memory structure was also
used. These methods ([3] and [6]) used existing methods to
compress the suffix arrays. Our method can be regarded as
one that boosts such methods in terms of memory storage and
query processing time.
2) Trajectory compression: As noted in Section I, shortestpath encodings have been used to compress spatial paths in
several papers [1], [2], [4] and [24]. As implied by Krumm
[29], an NCT dataset is expected to have a small k-th order
empirical entropy (Eq. (4)), but none of such shortest-pathbased compressors have provided an information-theoretic
evaluation of the compressed size. As an NCT compressor,

our method first focuses on high-order entropy and gives an
information-theoretic bound (Theorem 4). One of the methods
proposed in [1], MEL, is a different type of spatial path
compressor that achieves higher compressibility than shortestpath-based methods. As shown in Section V-D, RML achieves
a smaller entropy than MEL does. In [5], graph partitioning
was used to reduce the size of spatial paths.
To compress timestamps in NCTs, lossy compression methods are used in [1], [4], [5] and [24], whereas lossless
compression was used in [3]. These NCT compressors support
some useful queries. In any case, timestamps should be stored
and indexed to fit the purpose of use. In order to realize such
queries in a smaller storage, it is an interesting research direction to combine these timestamp compressors with CiNCT.
3) FM-index: FM-index, a compressed representation of
suffix arrays [7], was proposed by Ferragina and Manzini [8].
We have already described FM-index and the related topics in
Section II. Although there are a number of FM-index variants
(e.g., [16], [17], [20], [21]), these are essentially designed
for general strings. In Section VI, we compared our method
also with FM-indexes designed for large alphabet [20], [21].
For large σ, these methods can process suffix range queries
in O(|P | log log σ) time, which is much faster than typical
O(|P | log σ) time. Importantly, we employed the domainspecific knowledge of the target data (i.e., sparse transition
in road networks) to enhance the compression and query
processing. This point is the largest difference from the FMindex family designed for general strings.
VIII. C ONCLUSION
In this paper, we proposed CiNCT, a novel compressed
data structure capable of handling a very large number of
NCTs. We incorporated the sparsity of road networks into the
FM-index by using our proposed RML and PseudoRank techniques. The resulting data structure supports pattern matching
(i.e., via suffix range queries) and sub-path extraction from an
arbitrary position while still achieving high compressibility.
Our experiments showed that CiNCT outperformed existing
methods in terms of index size and search time, as shown
above Fig. 10, Table IV, and Fig. 12. Our method was even
faster than an uncompressed index. We also discussed theoretically why CiNCT is compact and fast. Further, we proved the
optimality of RML, i.e., the smallest size and the fastest search
are achieved. We also showed that RML performed even better
than the state-of-the-art NCT labeling method (MEL).
Our data structure has a wide range of applications in which
pattern matching based on spatial paths is a key component.
In fact, our method can be directly applied to some pioneering
methods for spatio-temporal NCT processing [3], [6]. Given
our method, we also expect that a practical spatio-temporal
database system will become possible in the future.

R EFERENCES
[1] Y. Han, W. Sun, and B. Zheng, “COMPRESS: A comprehensive
framework of trajectory compression in road networks,” ACM Trans.
Database Syst., vol. 42, no. 2, pp. 11:1–11:49, May 2017.
[2] G. Kellaris, N. Pelekis, and Y. Theodoridis, “Map-matched trajectory
compression,” J. Syst. Softw., vol. 86, no. 6, pp. 1566–1579, 2013.
[3] N. R. Brisaboa, A. Fariña, D. Galaktionov, and M. A. Rodrı́guez,
“Compact trip representation over networks,” in Proc. SPIRE’16, 2016,
pp. 240–253.
[4] B. Krogh, C. S. Jensen, and K. Torp, “Efficient in-memory indexing of
network-constrained trajectories,” in Proc. GIS ’16, no. 17, 2016.
[5] I. Sandu Popa, K. Zeitouni, V. Oria, and A. Kharrat, “Spatio-temporal
compression of trajectories in road networks,” GeoInformatica, vol. 19,
no. 1, pp. 117–145, 2014.
[6] S. Koide, Y. Tadokoro, and T. Yoshimura, “SNT-index: Spatio-temporal
index for vehicular trajectories on a road network based on substring
matching,” in Proc. SIGSPATIAL UrbanGIS Workshop’15, 2015.
[7] U. Manber and G. Myers, “Suffix arrays: a new method for on-line
string searches,” in Proc. SODA’90, 1990.
[8] P. Ferragina and G. Manzini, “Opportunistic data structures with applications,” in Proc. FOCS’00, 2000.
[9] M. Burrows and D. J. Wheeler, “A block-sorting lossless data compression algorithm,” in Technical Report 124.
[10] R. Grossi, A. Gupta, and J. S. Vitter, “High-order entropy-compressed
text indexes,” in Proc. SODA’03, 2003, pp. 841–850.
[11] G. Jacobson, “Space efficient static trees and graphs,” in Proc. FOCS’89,
1989, pp. 549–554.
[12] R. Raman, V. Raman, and S. Rao, “Succinct indexable dictionaries with
applications to encoding k-ary trees and multisets,” in Proc. SODA’02,
2002, pp. 233–242.
[13] G. Navarro, “Wavelet trees for all,” in Proc. CPM’12, 2012, pp. 2–26.
[14] V. Mäkinen and G. Navarro, “New search algorithm and time/space
tradeoffs for succinct suffix array,” in Tech.Rep. C-2004-20, Univ. of
Helsinki, 2004.
[15] G. Manzini, “An analysis of the Burrows-Wheeler transform,” J. ACM,
vol. 48, no. 3, pp. 407–430, 2001.
[16] J. Kärkkäinen and S. J. Puglisi, “Fixed block compression boosting in
FM-Indexes,” in Proc. SPIRE’12, 2011, pp. 174–184.
[17] V. Mäkinen and G. Navarro, “Implicit compression boosting with
applications to self-indexing,” in Proc. SPIRE’07, 2007, pp. 229–241.
[18] F. Claude and G. Navarro, “The wavelet matrix,” in Proc. SPIRE’12,
2012, pp. 167–179.
[19] G. Navarro and E. Providel, “Fast, small, simple rank / select on
bitmaps,” in Proc. SEA’12, 2012, pp. 295–306.
[20] A. Golynski, J. I. Munro, and S. S. Rao, “Rank/select operations on
large alphabets: A tool for text indexing,” in Proc. SODA’06, 2006, pp.
368–373.
[21] J. Barbay, T. Gagie, G. Navarro, and Y. Nekrich, “Alphabet partitioning
for compressed rank/select and applications,” in Proc. ISAAC’10, 2010,
pp. 315–326.
[22] S. Gog, A. Moffat, and M. Petri, “CSA++: fast pattern search for large
alphabets,” in Proc. ALENEX’17, 2017, pp. 73–82.
[23] N. Larsson and A. Moffat, “Offline dictionary-based compression,” in
Proc. DCC ’99, 1999, pp. 296–305.
[24] R. Song, W. Sun, B. Zheng, and Y. Zheng, “PRESS: A novel framework
of trajectory compression in road networks,” PVLDB, vol. 7, no. 9, pp.
661–672, 2014.
[25] P. Newson and J. Krumm, “Hidden Markov map matching through noise
and sparseness,” in Proc GIS’09, 2009, pp. 336–343.
[26] L. Nguyen-Dinh, W. G. Aref, and M. F. Mokbel, “Spatio-temporal access
methods: Part 2 (2003–2010),” IEEE Data Eng. Bull., vol. 33, no. 2, pp.
46–55, 2010.
[27] V. T. de Almeida and R. H. Güting, “Indexing the trajectories of moving
objects in networks,” Geoinformatica, vol. 9, no. 1, pp. 33–60, 2005.
[28] B. Krogh, N. Pelekis, Y. Theodoridis, and K. Torp, “Path-based queries
on trajectory data,” in Proc. GIS’14, 2014, pp. 341–350.
[29] J. Krumm, “A Markov model for driver turn prediction,” in Society of
Automotive Engineers (SAE) 2008 World Congress, 2008.
[30] “On compressing permutations and adaptive sorting,” Theoretical Computer Science, vol. 513, pp. 109–123, 2013.

A PPENDIX
A. Proof of Theorem 3
To begin with, let us introduce some mathematical notations.
Let us denote a set of integers as [σ] := {1, · · · , σ}. Consider
σ discrete probability distributions p1 , · · · , pσ on [σ] defined
by
nww0
,
(15)
pw0 (w) =
n·w0
where
nww0 is the number of bigrams ww0 in T and n·w0 =
P
0
w nww . First, we define a permutation of a distribution.
Definition 4: Let p be a discrete distribution on [σ]. A
permutated distribution pπ is a distribution where pπ (k) =
p(π(k)). Here π is a permutation on [σ].
In addition, we introduce the concept of a decreasing
distribution:
Definition 5: A discrete distribution p is decreasing iff
p(w) ≥ p(w + 1) for ∀w ∈ [σ]. Let F be a set of decreasing
distributions and F c be a set of non-decreasing distributions.
Note that we can always find a permutation π that makes any
distribution p decreasing, that is, pπ ∈ F.
Let us relate the above definitions to our problem. Since
any possible RML corresponds to an assignment of distinct
integers cww0 ∈ [σ] as mentioned in Section III-B, we can
regard it as an array of permutations Π = [π1 , · · · , πσ ]. We
denote such an encoder as φΠ . Our strategy, sorting by bigram
nww0 , corresponds to an array of permutations Π such that
each πi makes the distribution pi decreasing. Note that, if
w∈
/ Nout (w0 ), we can treat such cases as pw0 (w) = 0.
Our problem is to find an encoder φΠ that achieves the
minimum H0 (φΠ (Tbwt )). Consider a mixture distribution
X
pΠ =
αi pπi i
(16)
i∈[σ]

P

where αi = n·i / j n·j . Since the entropy of a discrete
distribution is defined as
X
H(p) = −
p(k) lg p(k),
(17)

Proof: Since g(x) = (x + ε) lg(x + ε) − x lg x is a
strictly increasing function, we have g(a) − g(b) > 0, which
is equivalent to Eq. 20.
Now we are ready to prove Theorem 7.
Proof of Theorem 7: We prove optimality by contradiction. Let Π+ be a set of permutations that minimizes H.
+
As discussed above, we can assume pΠ ∈ F without loss of
generality. Let us assume that there exists at least one πi ∈ Π+
such that pπi i ∈ F c . Let us define q := pπi i .
Since q ∈ F c , there exists k ∈ [σ] such that q(k) < q(k+1).
+
Based on Eq. 16, pΠ can be decomposed as
+

pΠ = (1 − αi )p̂ + αi q

(21)

P
πj
1
where p̂ := 1−α
j6=i αj pj . If p̂(k) ≤ p̂(k + 1), we have
i
+
+
pΠ (k) < pΠ (k + 1), which contradicts the assumption
+
pΠ ∈ F. Therefore, we have p̂(k) > p̂(k + 1).
Consider a permutation sk which only swaps k and k+1 and
a permutated distribution q sk . Let us define α := αi and β :=
1 − αi . We can calculate the difference of entropy functions
+
between the optimal solution pΠ = β p̂+αq and the swapped
sk
distribution β p̂ + αq :
+

H(pΠ ) − H(β p̂ + αq sk ) =
− {β p̂(k) + αq(k)} lg{β p̂(k) + αq(k)}
− {β p̂(k + 1) + αq(k + 1)} lg{β p̂(k + 1) + αq(k + 1)}
+ {β p̂(k + 1) + αq(k)} lg{β p̂(k + 1) + αq(k)}.

the following equality holds:
(18)

Therefore, we can reformulate our optimization problem as
follows.
Π∗ = argminΠ H(pΠ ).

−a lg a − (b + ε) lg(b + ε) + (a + ε) lg(a + ε) + b lg b > 0.
(20)

+ {β p̂(k) + αq(k + 1)} lg{β p̂(k) + αq(k + 1)}

k∈[σ]

H(pΠ ) = H0 (φΠ (Tbwt )).

Theorem 7: The optimal solution Π∗ consists of permutations such that each πi ∈ Π∗ makes the distribution pi
decreasing: pπi i ∈ F for ∀i ∈ [σ].
We first consider the following Lemma which qualitatively
implies that a more concentrated distribution has smaller
entropy.
Lemma 1: If a > b ≥ 0 and ε > 0, we have

(22)

Using the notation a = β p̂(k)+αq(k), b = β p̂(k +1)+αq(k),
and ε = αq(k + 1) − αq(k), we have a > b ≥ 0 and ε > 0.
Now Eq. 22 can be rewritten as
+

(19)

Consider an optimal Π∗ and any permutation π. Permutating
elements in Π∗ by π also yields another optimal solution
∗
∗
by definition: H(pΠ ) = H(pπ◦Π ) where π ◦ Π∗ = {π ◦
π1 , · · · , π ◦πn }. Here g ◦f indicates a composite function. We
∗
can therefore assume pΠ is a decreasing distribution without
loss of generality.
We now prove the following Theorem that directly leads to
Theorem 3.

H(pΠ ) − H(β p̂ + αq sk ) =
− a lg a − (b + ε) lg(b + ε) + (a + ε) lg(a + ε) + b lg b > 0.
where the last inequality holds from Lemma 1. This inequality
+
indicates that H(pΠ ) > H(β p̂ + αq sk ). That is,
[π1 , · · · , πi−1 , sk ◦ πi , πi+1 , · · · , πσ ]
is better than Π+ . However, this contradicts the optimality of
Π+ .

B. Proof of Theorem 4
We can prove Theorem 4 in a similar way to [17], which
proves the theorem for ICB with a balanced wavelet tree.
Proof: To begin with, we introduce some facts about
RRR [19]. Let us consider a bit vector B of length n.
The RRR divides B into small blocks of length b: B =
B (1) B (2) · · · B (n/b) . Each B (j) is represented by its class cj
and offset oj . Here the class is the number of 1’s in B (j) ,
and the offset is an index to distinguish the positions of 1’s in
B (j) . In fact, the total space needed for the classes becomes
the second term of Eq. 11 (see [19]). Since this term is already
considered in Eq. 12, what we have to evaluate is the
 offsets.
Each offset requires lg cbj bits because there are cbj possible
layouts of 1’s for the class cj .
Let us consider the partition of contexts of length k:
φ(Tbwt ) = L1 L2 · · · Ll (l ≤ σ k ). Since a bit vector in a node
v of a wavelet tree keeps the ordering, the bit vector Bv can
be divided into l blocks: Bv = B1v B2v · · · Blv . Here each Bjv
corresponds to Lj .
Now, we can consider small blocks of RRR which
is fully included in Bjv . Let us denote such blocks as
B (1) B (2) · · · B (t) . The offsets for these blocks require
 
t
X
b
lg
≤ |Bjv |H0 (Bjv )
(23)
c
j
j=1
bits. There are at most two blocks at the boundary of Bjv not
considered above. Their offset requires O(b) bits; therefore,
the offsets for Bjv need |Bjv |H0 (Bjv ) + O(b) bits in total.
Let us consider the space needed for Lj . Since there are at
most σ−1 inner nodes in a wavelet tree, summing the required
spaces over v, we obtain
X
|Bjv |H0 (Bjv ) + O(σb) = |Lj |H0 (Lj ) + O(σb) (24)
v

bits; the right hand side can be obtained by the recursive
calculation technique discussed in [30].
Summing the above equation over j, we have
l
X

|Lj |H0 (Lj ) + O(lσb).

(25)

j=1

Although Lj is an encoded string, the elements have a oneto-one correspondence with the non-encoded string because
of the definition of the RML. Hence, H0 (Lj ) is equal to
H0 (TW ), where W is the corresponding context and TW is
defined in Eq. 4.

