One of the basic problems in text searching is the *pattern matching problem*, which is also called the string matching problem. This is the problem of finding the occurrences of keywords (*patterns*) in a target text. Broadly speaking, there are two approaches to solving this problem. One is to build an index for the input text in advance, which is called *text indexing*. A text index allows for efficient searching when the target text is static and is not subsequently updated. The other is to access the input text sequentially from the beginning to the end while checking if the given pattern matches at the current reference position in the text. This is called *text scanning*. Text scanning is applicable even if the text is updated from time to time and it does not require index structures. In general, text indexing is superior in terms of search speed, while text scanning is superior in terms of search flexibility. By convention, "pattern matching" often refers to the latter, text scanning, in a narrow sense.

In this chapter, we outline a fusion technology of data compression and pattern matching called *compressed pattern matching*. First, in Sect. 7.2, we look back on the history of this field of study. Then, in Sect. 7.3, we provide some notation and definitions that are used in the following sections. In addition, we recall grammar compression, which is the key compression scheme for compressed pattern matching. Next, in Sect. 7.4, we introduce the general framework of compressed pattern matching proposed by Kida et al. [21]. In Sects. 7.5 and 7.6, we present outlines of Repair-VF and MR-Repair, respectively, which are the results of our work in this study. Finally, we conclude in Sect. 7.7.

#### **7.2 History of Compressed Pattern Matching Research**

The technology of combining data compression and pattern matching emerged in the early 1990s. This has come to be known as the *compressed pattern matching problem* [1], which is the problem of performing pattern matching without first decompressing the compressed input text. Formally, when a text *T* = *t*1*t*2 ... *tu* (*ti* is a symbol) is given in compressed form *Z* = *z*1*z*2 ...*zn* (*zi* is an element of the compressed text), and pattern *P* is given, the problem is to find all occurrences of *P* in *T* using only *Z* and *P*. A simple method is to first decompress *Z* to *T* and then use some commonly used pattern matching algorithm. However, this approach requires O(*m* + *u*) time for pattern matching in addition to the decompression time.

The optimal algorithm for the compressed pattern matching problem is one that performs pattern matching in O(*m* + *n*) time in the worst case. However, it is not easy to achieve both efficient compression of text data and fast pattern matching on it. In the initial research in this field, individual pattern matching algorithms were developed for each compression method. For example, Eilam-Tzoreff and Vishkin et al. [15] proposed an algorithm for run-length compression, G¸asieniec et al. [18] and Farach Thorup et al. [16] proposed algorithms for LZ77, and Amir et al. [2] proposed an algorithm for LZW [54]. However, these algorithms tend to be complicated, and as Manber [29] pointed out, it is questionable as to whether they are more practical than the simple method.

From the late 1990s to the early 2000s, several efficient methods for compressed pattern matching emerged [22, 38, 41]. For the first time, it was shown experimentally that these methods can perform pattern matching faster than the simple method. Furthermore, methods have appeared that can perform pattern matching faster by about the compression ratio than matching on the original text. The key is to select a compression method suitable for pattern matching even at the expense of compression ratio. In fact, Byte Pair Encoding (BPE), which was used by Shibata et al. [48] for this purpose, has a compression ratio of at most about 50% for natural language texts, while the LZ-family methods can compress the same texts to about 30% or less. However, text compression by BPE offers an advantage for pattern patching because all the codewords are fixed at 8 bits and the correspondence between each codeword and a portion of the text is relatively clear.

This caused a paradigm shift. Whereas individual pattern matching algorithms were previously developed for each data compression method, we realized that in order to increase the matching speed it would be better to develop a new data compression method suitable for pattern matching. In fact, in the 2000s, several data compression methods were proposed for this purpose.

One of the main groups of compression methods based on this idea are the compression methods proposed by Brisaboa et al. [7–9] and by Klein and Ben-Nissan [24]. These are based on a technique called *dense coding* [8]. Dense coding divides an input (natural language) text into words, and then encodes them so that the codewords become shorter in descending order of the frequency of the words. In addition, each codeword is assigned a bit pattern that has an explicit end to facilitate codeword extraction. Although dense coding offers good performance in terms of both compression ratio and pattern matching speed, some ingenuity is required to apply it to data such as DNA sequences that cannot be divided into words.

The other system is *grammar-based compression* (or *grammar compression*) [23] with fixed-length coding. This idea is an extension of BPE and can be applied even if an input text cannot be separated into words. One direct improvement of BPE is a method using a context-sensitive grammar by Maruyama et al. [34], while for compression methods based on context-free grammar we have the methods by Klein and Shapira [25] and Uemura et al. [51]. In both methods, a modified version of suffix tree [53] is used as a dictionary tree for constructing grammar.

In this chapter, for convenience, we refer to the former system as the dense coding system and the latter as the VF coding system.

In the 2010s, new data compression algorithms began to appear that achieved compression performance comparable to well-known compression tools such as Gzip and Bzip while maintaining properties suitable for pattern matching. Among the two systems described above, the VF coding system has difficulties in terms of compression rate and compression speed. Therefore, research looked into searching for and improving grammar compression, which is the basis of the VF coding system. Among this work, the algorithm RePair [26], which was proposed before the name "grammar compression" was used, has attracted attention because of its excellent compression ratio. Yoshida and Kida et al. [56] proposes a variant of RePair, called Repair-VF, which reduces the decrease in compression ratio by suppressing unnecessary grammar rules while encoding the output using fixed length codewords. The time and space complexities required for Repair-VF are both O(*n*) for text of length *n*, which is the same as the original RePair. Repair-VF realizes high-speed pattern matching on compressed text while having a good compression ratio comparable to Gzip.

Very recently, we proposed a novel variant of RePair, called MR-RePair [17], which constructs more compact grammars than RePair, particularly for highly repetitive texts. This achievement comes from an analysis of RePair. We show in [17] that the main process of RePair, that is, the step by step substitution of the most frequent symbol pairs, works within the corresponding most frequent *maximal repeats*. We then reveal in [17] the relationship between maximal repeats and grammars constructed by RePair.

#### **7.3 Preliminaries**

#### *7.3.1 Definitions of Notation and Terms*

Let be an *alphabet*, that is, an ordered finite set of symbols. An element *T* = *t*1 ... *tn* of -∗ is called a *string* or a *text*, where |*T* | = *n* denotes its length. Let ε be an empty string of length 0, that is, |ε| = 0. We denote a concatenation of two strings *x*, *y* ∈ - by *x* · *y*, or *x y* for simplicity if no confusion occurs.

If *T* = *xyz* with *x*, *y*,*z* ∈ -∗, then *x*, *y*,*z* are called a *prefix*, *substring*, and *suffix* of *T* , respectively. Let *T* [*i* : *j*] = *ti* ··· *tj* for any 1 ≤ *i* ≤ *j* ≤ *n* denote a substring of *T* beginning at *i* and ending at *j* in *T* , and let *T* [*i*] = *ti* denote the *i* th symbol of *T* . Let w[*i* : *j*] = ε if *j* < *i* for simplicity.

#### *7.3.2 Grammar Compression*

A *context-free grammar* (*CFG* or simply *grammar*) G is defined as a four-tuple *G* = {*V*,-, *S*, *R*}, where *V* denotes an ordered finite set of variables, denotes an ordered finite alphabet, *R* denotes a finite set of binary relations called *production rules* (or *rules*) between *V* and (*V* ∪ -)∗, and *S* ∈ *V* denotes a special variable called *start variable*. A production rule refers to the situation where a variable is substituted and written in the form v → w, with v ∈ *V* and w ∈ (*V* ∪ -)∗. Let *X*, *Y* ∈ (*V* ∪ -)∗. If there are *xl*, *x*, *xr*, *y* ∈ (*V* ∪ -)∗ such that *X* = *xl x xr*, *Y* = *xl yxr*, and *x* → *y* ∈ *R*, we write *X* → *Y* , and denote the reflexive transitive closure of→as ∗ ⇒. Let *val*(v) be a string derived from v, that is, v ∗ ⇒ *val*(v). We define grammar *G*ˆ = {*V*ˆ , -, ˆ *S*ˆ, *R*ˆ} as a *subgrammar* of *G* if *V*ˆ ⊆ *V*, -ˆ ⊆ (*V* ∪ -), and *R*ˆ ⊆ *R*.

Given a text *T* , *grammar compression* is a method for lossless text data compression that constructs a restricted CFG uniquely deriving the text *T* . For *G* to be deterministic, the production rule for each variable v ∈ *V* must be unique. In what follows, we assume that every grammar is deterministic and each production rule is v*i* → *expri* , where *expri* is an expression either *expri* = *a* (*a* ∈ -) or *expri* = v*j*1 v*j*2 ...v*jn* (*i* > *jk* for all 1 ≤ *k* ≤ *jn*). To estimate the effectiveness for compression, we use the size of the constructed grammar, which is defined as the total length of the right-hand side of all production rules of the grammar.

While the problem of constructing the smallest such grammar for a given text is known to be NP-hard [10], several approximation algorithms have been proposed. One of them is RePair [26], which is an off-line grammar compression algorithm. Despite its simple scheme, RePair is known for its high compression in practice [12, 19, 52], and hence, it has been comprehensively studied. Some examples of studies on the RePair algorithm include its extension to an online algorithm [35], practical working time/space improvements [6, 47], applications to various fields [12, 27, 49], and theoretical analysis of generated grammar sizes [10, 40, 42].

#### **7.4 Framework for Compressed Pattern Matching**

A grammar compressed text can be expressed in a framework called collage systems [21]. A pattern matching algorithm on the compressed text can then be obtained as an instance of the general algorithm on the collage system. Algorithm on collage systems can be understood as an extension of the Knuth-Morris-Pratt method (KMP method) [14].

#### *7.4.1 KMP Method*

The KMP method a well-known linear-time algorithm for pattern matching on an ordinary (uncompressed) text. Its movement can be modeled as a linear automaton (KMP automaton) for a given pattern *P*.

For a given pattern *P*, a KMP automaton consists of two functions:

goto function*g* : *Q* × - → *Q* ∪ {*fail*}, failure function *f* : *Q* \ {0} → *Q*,

where *Q* = {0, 1,..., |*P*|} is the set of states, and *fail* is a special symbol that is not included in *Q*. For *j* ∈ *Q* and *a* ∈ -, the goto function *g* returns *j* + 1 if *P*[*j* + 1] = *a* holds, otherwise it returns*fail*. For *j* = 0, let *g*(0, *a*) = 0 for all *a* ∈ - where *P*[1] = *a* holds. For *j* ∈ *Q* \ {0}, the failure function *f* returns the maximum integer *k* such that *P*[1 : *k*] = *P*[*j* − *k* + 1 : *j*] holds.

**Fig. 7.1** KMP automaton for *P* = abacb. In this figure, each circle indicates a state, and the double circle indicates the final state. Solid arrows and dashed arrows indicate the goto function and failure function, respectively

**Fig. 7.2** Movement of KMP automaton. Solid arrows and dashed arrows indicate state transitions caused by the goto function and the failure function, respectively

The automaton repeats state transitions by tracing *g* corresponding to the characters read one by one from the input text. If *g* returns *fail*, then *f* is repeatedly called with the current state number to go back until a transition by *g* succeeds with the same character. When the automaton finally reaches the rightmost state, it can be judged that *P* has occurred.

Figure 7.1 shows the KMP automaton for pattern *P* = abacb. The movement of the KMP automaton in Fig. 7.1 for the text *T* = abacbbabaabacb is shown in Fig. 7.2. In this example, it can be judged that *P* has occurred when the automaton reaches the state number 5.

To eliminate the failure function, we define the state transition function δ : *Q* × -→ *Q* as follows:

$$\delta(j, a) = \begin{cases} g(j, a) & \text{if } g(j, a) \neq fail, \\ \delta(f(j), a) \text{ otherwise} \end{cases}$$

Moreover, we extend it to *Q* × -∗ as follows:

$$
\delta^\*(j,\varepsilon) = j, \quad \delta^\*(j,\mu a) = \delta(\delta^\*(j,\mu), a),
$$

where *j* ∈ *Q*, *u* ∈ -∗, and *a* ∈ -.

#### *7.4.2 Collage System*

A *collage system* is a pair D,S defined as follows: D is a sequence of assignments *X*1 = *expr*1; *X*2 = *expr*2;··· ; *Xn* = *exprn*, where each *Xk* is a token and *exprk* is any of the form:

*a* for *a* ∈ - ∪ {ε}, (*primitive assignment*) *Xi X j* for *i*, *j* < *k*, (*concatenation*) [*j*] *Xi* for *i* < *k* and an integer *j*, (*prefix truncation*) *X*[*j*] *i* for *i* < *k* and an integer *j*, (*suffix truncation*) (*Xi*)*j* for *i* < *k* and an integer *j*. (*j times repetition*)

Let the set of all tokens in D be denoted by *F*(D). Each token represents a string obtained by evaluating the expression as it implies. Let the string represented by token *X* ∈ *F*(D) be denoted by *X*.*u*. For example, for *X*1 = a; *X*2 = b; *X*3 = *X*1 · *X*2; *X*4 = (*X*3)3; *X*5 = *X*[1] 4 , *X*4.*u* = ababab and *X*5.*u* = ababa. However, in this section we identify token *X* with the string it represents, and simply denote both by *X* unless confusion occurs.

Let the number of assignments in D be the *size* of D, and denote it by ||D||, that is, ||D|| = |*F*(D)| = *n*. For a sequence S = *Xi*1 , *Xi*2 ,..., *Xik* of tokens defined in D, we denote by |S| the number *k* of tokens in S.

The collage system D,S represents the string obtained by concatenating *Xi*1 ,..., *Xik* . That is, D and S correspond to the dictionary and compressed text in a compression method, respectively. Both D and S can be encoded in various ways. The compression ratios therefore depend on their encoding sizes rather than ||D|| and |S|.

A collage system is said to be *truncation-free* if D contains no truncation operation. A collage system is said to be *regular* if D contains neither truncation nor repetition operations. A regular collage system is said to be *simple* if for every assignment *X* = *Y Z*, |*Y*.*u*| = 1 or |*Z*.*u*| = 1.

#### *7.4.3 Pattern Matching on Collage Systems*

The basic idea of pattern matching on a collage system is to simulate the movement of the KMP automaton on uncompressed text. Using the state transition function δ∗ of the KMP automaton defined in Sec. 7.4.1, we define the function *Jump* : *Q* × *F*(D) → *Q* as follows:

$$
\
Mump(j, X) = \delta^\*(j, X.u).
$$

The intent of *Jump* is to simulate the state transition of the original KMP automaton by jumping when it receives token *X*. Moreover, for any *j* ∈ *Q* and *X* ∈ *F*(D), we define the set *Output*(*j*, *X*) = *OccP* (*P*[1 : *j*] · *X*.*u*), where *OccP* (*x*) indicates the set of all indices of occurrences of *P* within *x*.

**Fig. 7.3** A matching algorithm on a collage system

For a given collage system D,S representing text *T* and for a given pattern *P*, the pattern matching algorithm preprocesses the information required to calculate *Jump* and *Output* from D, and then performs matching while scanning a sequence of tokens in S one by one from the head (Fig. 7.3).

From the results of [21], the following theorem is obtained.

**Theorem 1** (Theorem 3 of [21]) *For a collage system* D,S*, the compressed pattern matching problem can be solved in* O(||D|| + |S| + *m*2 + *R*) *time using* O(||D|| + *m*2) *space if* D,S *is regular, where R is the number of occurrences of pattern P in the text represented by* D,S*.*

This theorem applies to both RePair and Repair-VF because texts compressed by these can be described by regular collage systems.

#### **7.5 Repair-VF**

This section first introduces RePair and then gives an outline of Repair-VF. Repair-VF has a structure that combines RePair with a fixed-length coding. Please refer to the literature [56] for the details of Repair-VF and experimental results for its performance.

#### *7.5.1 RePair*

RePair is a grammar compression algorithm that was proposed by Larsson and Moffat [26]. For input text *T* , let *G* = {*V*,-, *S*, *R*} be the grammar constructed by RePair. The RePair procedure can then be described by the following steps:

**Step 1**. Replace each symbol *a* ∈ with a new variable v*a* and add v*a* → *a* to *R*. **Step 2**. Find the most frequent pair *p* in *T* .

**Fig. 7.4** Example of the grammar generation process of RePair for *T* = abracadabra. The generated grammar is {{va, vb, vr, vc, v1, v2, v3, *S*},{a,b,r,c,d}, *S*,{va → a, vb → b, vr → r, vc → c, vd → d, v1 → vavb, v2 → v1vc, v3 → v2vd, *S* → v3vcvavdv3}} with a size of 16

Figure 7.4 illustrates an example of the grammar generation process of RePair.

The following theorem relates to the performance of RePair shown by Larsson and Moffat [26].

**Theorem 2** ([26]) *RePair works in* O(*n*) *expected time and* 5*n* + 4*k*2 + 4*k* + √*n* + 1 − 1 *words of space, where n is the length of the source text, k denotes the cardinality of the source alphabet, and k denotes the cardinality of the final dictionary.*

#### *7.5.2 Outline of Repair-VF*

The original RePair encodes the rules in *R* excluding *S* using Elias gamma coding, that is, each codeword has a variable length, whereas Repair-VF uses a fixed-length code. The right side of *S* corresponds to the compressed text, and is converted to a sequence of fixed-length codewords of the rules in *S*.

Consider the number of fixed-length coded rules. In the process of Step 1 of RePair, |-| rules are created. In addition, the process of Step 3 of RePair replaces the most frequent pair and at the same time adds one rule to *R*. Let *s* be the number of rules which are added to *R* in Step 3. The total number of rules is then |-| + *s*, and thus each symbol can be fixed-length encoded with log(|-| + *s*) bits. The information about can be restored from the rules added in Step 3, so there is no need to explicitly save it. Therefore, only the rules added in Step 3 and the right side of *S* added last in Step 4 need to be saved. In the former, the right side of each rule consists of two symbols, so the total number of symbols to be saved is 2*s*. Since the latter depends on the input text *T* , let *n* be the length of *T* . The number of bits of the output compressed data is then

$$(2s+n)\lceil \log(|\Sigma|+s) \rceil\tag{7.1}$$

bits in total.

We want to find the best*s* that minimizes the total number of output bits. Note that the final output tends to be smaller as*s* increases up to some point. In RePair, Step 3 is repeated until the frequency of the most frequent pair becomes 1. In the case of using a fixed-length code as above, this increases useless rules. Increasing the number of rules can increase the bit length per symbol, resulting in a longer final bit length. We can eliminate the waste by terminating the process in the middle. However, it is difficult to determine on the fly the *s* at which the output becomes minimal. Even after the first time the output size increases, the length of *S* may become shorter by continuing Step 3, and the output size may decrease again.

Therefore, during the processing of RePair, we record the minimum value of the output size and the corresponding *s* by calculating Equation (7.1) every time a rule is added in Step 3. Note that the calculation of Equation (7.1) does not require actual coding or outputting since a fixed-length code is used. Finally, when *S* is output in Step 4, we can obtain the smallest output by outputting while expanding the rules added after the best *s*.

This is *Repair-VF* (called Repair-VF-best in the original paper). The suffix "VF" comes from an abbreviation for variable-to-fixed length coding (VF coding). For the input text *T* , each rule of the output grammar *G* corresponds to a substring of *T* , and the right-hand side of *S* can be regarded as the variable length factorization of *T* . Thus, Repair-VF can be viewed as a VF coding from the viewpoint of information source coding.

#### **7.6 MR-Repair**

In this section we outline MR-Repair, which is a method to reduce the output grammar size by focusing on the relationship between RePair and maximal repeats. Please refer to the literature [17] for the details of MR-Repair and experimental results for its performance.

#### *7.6.1 Maximal Repeats*

Let *s* be a substring of text *T* . If the frequency of *s* is greater than 1, *s* is called a *repeat*. A *left* (or *right*) *extension* of *s* is any substring of *T* in the form w*s* (or *s*w), where w ∈ -∗. We define *s* as a *left* (or *right*) *maximal* if left (or right) extensions of *s* occur a strictly lower number of times in *T* than *s*. Accordingly, *s* is a maximal repeat of *T* if *s* is both left and right maximal. In this paper, we only consider strings with a length of more than 1 as maximal repeats. For example, the substring abra of *T* = abracadabra is a maximal repeat, whereas br is not.

The following theorem describes an essential property of RePair, that is, RePair recursively replaces the most frequent maximal repeats.

**Theorem 3** (Theorem 1 of [17]) *Let T be a given text, under the condition that every most frequent maximal repeat of T does not appear overlapping itself. Let f be the frequency of the most frequent pairs of T , and t be a text obtained after all pairs with frequency f in T are replaced by variables. There is then a text s such that s is obtained after all maximal repeats with frequency f in T are replaced by variables, and s and t are isomorphic to each other.*

#### *7.6.2 MR Order*

According to Theorem 1 of [17], if there is just one most frequent maximal repeat in the current text, then RePair replaces all occurrences of it step by step. However, a problem arises if there are two or more most frequent maximal repeats, with some of them overlapping. In this case, the selection order of pairs (of course, they are most frequent) affects the priority of maximal repeats. We call this order of selecting (summarizing) maximal repeats the *maximal repeat selection order* (or simply *MRorder*). Note that the selection order of pairs actually depends on the implementation of RePair. If there are several distinct most frequent pairs with overlaps, RePair constructs grammars with different sizes according to the selection order of the pairs.

However, the following theorem states that the MR-order rather than the replacement order of pairs determines the size of the grammar generated by RePair.

**Theorem 4** (Theorem 2 of [17]) *The sizes of grammars generated by RePair are the same if they are generated in the same MR-order.*

#### *7.6.3 Algorithm*

The main strategy of the proposed method is to recursively replace the most frequent maximal repeats instead of the most frequent pairs.

**Definition 1** (*Definition 3 of* [17]) For an input text *T* , let *G* = {*V*,-, *S*, *R*} be the grammar generated by MR-Repair. MR-Repair constructs *G* through the following steps:

**Step 1.** Replace each symbol *a* ∈ with a new variable v a and add v*a* → *a* to *R*.

**Step 2.** Find the most frequent maximal repeat *r* in *T* .

**Step 3.** Check if |*r*| > 2 and *r*[1] = *r*[|*r*|], and if so, use *r*[1 : |*r*| − 1] instead of *r* in Step 4.

**Fig. 7.5** Example of the grammar generation process of MR-Repair for *T* = abracadabra. The generated grammar is {{v*a* , v*b*, v*r*, v*c*, v*d* , v1, *S*},{*a*, *b*,*r*, *c*, *d*}, *S*,{v*a* → *a*, v*b* → *b*, v*r* → *r*, v*c* → *c*, v*d* → *d*, v1 → v*a* v*b*v*r*, v2 → v1v*a* , *S* → v2v*c*v*a* v*d* v2}} with a size of 15

Figure 7.5 shows an example of the grammar generation process of MR-Repair. As shown in this figure, the size of the grammar generated by MR-Repair is smaller than that generated by RePair shown in Figure reffig:repair.

**Theorem 5** (Theorem 5 of [17]) *Assume that RePair and MR-Repair work based on the same MR-order for a given text. Let grp and gmr be the sizes of the grammars generated by RePair and MR-Repair, respectively. Then,* 1 2 *grp* < *gmr* ≤ *grp holds.*

#### **7.7 Conclusion**

In this chapter, we outlined research on compressed pattern matching and showed that we can speed up pattern matching by selecting a suitable compression method. This has led to the development of compression methods that are useful for pattern matching. Whereas the initially developed compression methods had low compression performance, Repair-VF [56] achieves both a good compression rate and good matching speed by combining advanced grammar compression with fixed-length code. Collage systems [21] provide a unified algorithm for compressed pattern matching, allowing us to obtain an efficient pattern matching algorithm for grammar compression as an instance of the unified algorithm.

Since proposing Repair-VF, we have proposed several improvements for it. LT-Repair [35] improves RePair processing semi-online by adding the constraint called the *left-tall condition* to its grammar. This makes it possible to efficiently compress large-scale text data with small memory.

MR-Repair [17], which we have recently proposed, is a method that reduces the output grammar size by focusing on the relationship between RePair and maximal repeats. Although heuristic improvements [4, 20, 28, 36] focusing on non-maximal repetitive substrings have previously been proposed, MR-Repair is superior because it has been proven to generate theoretically smaller grammar than the original RePair. A topic for future work is to see whether compressed pattern matching using these methods can be performed efficiently.

In the present work, we mainly explained the compressed pattern matching problem based on text scanning. However, the *succinct index* technology which combines text index and data compression was also established in 2000. This is an indexing technology that utilizes a succinct data structure that can solve query processing with a small space of almost the information-theoretic lower bound. Succinct index has an excellent property that allows full-text searching while compressing a target text smaller than the original text. For details of this technology, refer to the excellent book by Navarro [37].

In terms of online grammar compression methods, there exists FOLCA proposed by Maruyama et al. [33] and its improvement SOLCA proposed by Takabatake et al. [50]. FOLCA is based on a string factorization called edit-sensitive parsing. It performs factorization and grammar generation in parallel while reading an input text sequentially from the beginning. It is known that *straight line programs* (restricted CFGs) generated by FOLCA can be used as index structures [5].

In recent years, Martinez et al. [31] proposed a novel compression method called Marlin and an improvement of it [32]. These methods achieve both decompression at ultrahigh-speed and good performance in terms of compression ratio. If we can decompress compressed data at sufficiently high speed, we can perform pattern matching efficiently even if it is performed after decompressing the data. Comparative studies on these approaches are also left for future work.

#### **References**

**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

# **Chapter 8 Orthogonal Range Search Data Structures**

**Kazuki Ishiyama and Kunihiko Sadakane**

**Abstract** We first review existing space-efficient data structures for the orthogonal range search problem. Then, we propose two improved data structures, the first of which has better query time complexity than the existing structures and the second of which has better space complexity that matches the information-theoretic lower bound.

#### **8.1 Introduction**

Consider a set *P* of *n* points in the *d*-dimensional space R*d* . Given an orthogonal range *Q* = - *l* (*Q*) 0 , *u*(*Q*) 0 × - *l* (*Q*) 1 , *u*(*Q*) 1 ×···× - *l* (*Q*) *d*−1, *u*(*Q*) *d*−1 , the problem of answering queries for information on *P* ∩ *Q*, the subset of *P* contained in the range *Q*, is called the *orthogonal range search* problem, and is one of the fundamental problems in computational geometry.

The information obtained about *P* ∩ *Q* differs depending on the query. The most basic queries are the *reporting* query, which enumerates all the points in *P* ∩ *Q*, and the *counting* query, which returns the number of points |*P* ∩ *Q*|. There are other queries such as the *emptiness* query, which checks whether *P* ∩ *Q* is empty or not, and aggregate queries, which compute the summation, average, or variance of weights of points in the query range.

Applications of the orthogonal range search problem include database searches [21]. For example, assuming there is a database of employees of a company, then a query to count the number of employees whose duration of service is at least *x*1 years and at most *x*2 years, age is at least *y*1 and at most *y*2, and annual income is at least *z*1 and at most *z*2, can be formalized as an orthogonal range search problem. Other applications include geographical information systems, CAD, and computer graphics.

K. Ishiyama · K. Sadakane (B)

Graduate School of Information Science and Technology, The University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-8656, Japan e-mail<EMAIL_ADDRESS>
© The Author(s) 2022

N. Katoh et al. (eds.), *Sublinear Computation Paradigm*, https://doi.org/10.1007/978-981-16-4095-7\_8

In such applications, it is common to perform multiple queries on the same point set *P*. We therefore consider constructing the problem as an indexing problem: Given a point set *P* a priori, we first construct some data structure *D* from *P*. Then, when a query range *Q* is given, we answer the query using the data structure *D*.

#### *8.1.1 Existing Work*

In many existing works, the number *n* of points is regarded as a variable for evaluating time complexity and the number *d* of dimensions is regarded as a constant. However, in this chapter, we regard *d* as a variable too. For the computation model, we use w-bit word RAM where w = - (lg *n*) bits. That is, a constant number of coordinate values can be treated in constant time. Then, it takes O(*d*) time to check whether a point is inside a query range.

If more space than - (*dn*) words is allowed to be used for the space complexity of data structures and if we assume that *d* is a constant, then we can perform the counting and reporting queries in time polynomial to log *n*. *Range trees* [2, 14, 15, 23] are such data structures. Range trees support counting queries in O *d* log*d*−1 *n* time and reporting queries in O *d* log*d*−1 *n* + *dk* time using O *dn* log*d*−1 *n* word space, where *k* = |*P* ∩ *Q*|, that is, the number of points enumerated by a reporting query using the fractional cascading technique [15, 23]. Although these data structures are time-efficient, it is desirable to develop more space-efficient data structures.

Some data structures having linear space complexity have been proposed. For example, quad trees [6] were the first data structures used for orthogonal range search. Unfortunately, quad trees have terrible worst-case behaviors. To overcome this, kd-tree [1] is used. The query time complexity of the kd-tree is O *d*2*n d*−1 *d* for counting and O *d*2*n d*−1 *d* + *dk* for reporting [13].

These data structures store the coordinates of points separately in plain form, and therefore can be applied to the case of real-valued coordinates. However, if the coordinates take integer values from 0 to *n* − 1, then there exist data structures with even smaller space complexity and query time complexity. For example, Chazelle [4] proposed a data structure for the two-dimensional case with linear space complexity and time complexity of O(lg *n*)for counting and O(lg *n* + *k* lgε *n*)for reporting where 0 <ε< 1 is any constant. Note that although the assumption that each coordinate value is an integer from 0 to *n* − 1 seems too strict, as is explained in Sect. 8.2.2, any orthogonal range search problem in *d*-dimensional space can be reduced into one on the [*n*] *d* grid, and therefore the assumption does not create any difficulties.

There has also been research on succinct data structures for the orthogonal range search problem. The wavelet tree [9] is a data structure which was originally proposed for representing compressed suffix arrays, and it later turned out that wavelet tree can support various queries efficiently [18]. For the orthogonal range search problem, wavelet tree can support counting queries in O(lg *n*) time and reporting queries in O((1 + *k*)lg *n*) time [8]. Bose et al. [3] proposed improved succinct data structures that support counting queries in O(lg *n*/ lg lg *n*) time and reporting queries in O(((1 + *k*)lg *n*/ lg lg *n*) for two-dimensional cases.

For higher dimensions, Okajima and Maruyama [20] proposed the KDW-tree, which is a succinct data structure for any dimensionality. The query time complexity of the KDW-tree is smaller than that of the kd-tree. If we assume *d* is a constant, counting queries take O *n d*−2 *d* lg *n* time and reporting queries take O*n d*−2 *d* + *k* lg *n* time. The KDW-tree has been shown to be practical by numerical experiments.

#### *8.1.2 Our Results*

We show space and time complexities of data structures for the orthogonal range search problem explained in Sect. 8.1.1 and our proposed data structures in Table 8.1. Note that these are for the case where the coordinates are integers from 0 to *n* − 1, and the space complexities are measured in bits. Table 8.1 shows reporting time complexities. Counting time complexities can be obtained by letting *k* = 0.

Our data structures are space-efficient for high-dimensional orthogonal range search problems.

Our first data structure has the same space complexity as the KDW-tree and better query time complexities. Note that the result in Table 8.1 is for the case of *d* ≥ 3. If *d* = 2, we can improve the *n d*−2 *d* term to lg *n*. This result appeared in [11].

Note that, as shown in Sect. 8.2.1, the necessary space to represent a set of *n* points in *d*-dimensional space such that each coordinate takes an integer value from 0 to *n* − 1 is (*d* − 1)*n* lg *n* + - (*n*) bits. This means that if we assume *d* is a constant, the space complexity of the KDW-tree and our first data structure does not match the information-theoretic lower bound asymptotically.

**Table 8.1** Comparison of complexities. The results or KDW-tree and Ours 1 are for *d* ≥ 3. Note that *k* is the number of points enumerated by a reporting query. The time complexities for counting queries are obtained by letting *k* = 0 in the time complexities for reporting queries

Our second data structure uses (*d* − 1)*n* lg *n* + (*d* − 1) · o(*n* lg *n*) bits of space. This asymptotically matches the information-theoretic lower bound even if *d* is assumed to be a constant. Therefore, we can say this data structure is truly succinct. Unfortunately, the worst-case query time complexity is O(*dn* lg *n*), which is not fast in theory. However, this data structure is fast in practice for the case where the number *d* of dimensions is large but the number *d* of dimensions used for a query is small. This kind of query often occurs in the database search applications shown in Sect. 8.1. This result appeared in [10].

#### **8.2 Preliminaries**

In this paper, we assume that coordinates of points are non-negative integers. As will be explained in Sect. 8.2.2, we sometimes assume that coordinates are integers from 0 to *n* − 1. Therefore, we define [n] as the set {0, 1,..., *n* − 1}. For a d-dimensional space, we denote each dimension by dim. 0, dim. 1, . . . , dim. *d* − 1, coordinate values of a point by 0-th coordinate value, 1-th coordinate value, ..., *d* − 1-th coordinate value. For a rooted tree, we assume the depth of the root node is 0. Throughout the paper, log x denotes the natural logarithm and lg x denotes the base 2 logarithm.

Next, we define two concepts used in this chapter. The first one is containment degree. This is the concept of an inclusion relationship between two orthogonal ranges introduced in [20]. For two *d*-dimensional orthogonal ranges *Q* = - *l* (*Q*) 0 , *u*(*Q*) 0 × ... × - *l* (*Q*) *d*−1, *u*(*Q*) *d*−1 and *R*= - *l* (*R*) 0 , *u*(*R*) 0 ×···× - *l* (*R*) *d*−1, *u*(*R*) *d*−1 , we define CDeg(*R*, *Q*) as

$$\text{CDeg}(R,\,\mathcal{Q}) = \# \left\{ i \in [d] \; \Big|\; \left[l\_i^{(R)}, u\_i^{(R)}\right] \subseteq \left[l\_i^{(\mathcal{Q})}, u\_i^{(\mathcal{Q})}\right] \right\}$$

and call it the containment degree of *R* with respect to *Q*. This is the number of dimensions, in each of which *R* is contained in *Q*. The containment degree is an important concept for analyzing time complexities of orthogonal range search algorithms.

Next, we explain *z*-value. This is a projection of multi-dimensional data onto onedimensional data as proposed by Morton [17]. Consider a point *p* = (*p*0, *p*1,..., *pd*−1) in the *d*-dimensional space where the coordinate values are integers. If coordinate values are expressed as *l*-bit binary numbers *p*0 = *b*0 0*b*1 0 ··· *bl*−1 0 , *p*1 = *b*0 1*b*1 1 ··· *bl*−1 1 ,..., *pd*−1 = *b*0 *d*−1*b*1 *d*−1 ··· *bl*−1 *d*−1, the *z*-value *z*(*p*) of point *p* is defined as

$$z(p) = b\_0^0 b\_1^0 \cdots b\_{d-1}^0 b\_0^1 b\_1^1 \cdots b\_{d-1}^1 \cdots b\_0^{l-1} b\_1^{l-1} \cdots b\_{d-1}^{l-1} \cdots$$

In the case of a two-dimensional space, if we arrange grid points in increasing order of *z*-value, we see a z-shape curve as shown in Fig. 8.1. We therefore call the value *z*-value.

**Fig. 8.1** Curve obtained by joining grid points in *z*-valueorder in two-dimensional space

## *8.2.1 Succinct Data Structures and Information-Theoretic Lower Bound*

Succinctness of data structures was proposed by Jacobson [12] and is one of the criteria for measuring space complexities of data structures. It is defined as follows.

Let *n* be the number of different values that an object can take. Then, we need at least lg *n* bits of space to represent the object. If the space complexity *S*(*n*) of a data structure representing the object satisfies *S*(*n*) = lg *n* + o(lg *n*) bits, we say the data structure is *succinct* and lg *n* bits is the *information-theoretic lower bound* of the size of representations of the object. Note that succinct data structures not only offer data compression, but also support some efficient queries. For orthogonal range search, a naive algorithm supports linear time queries by scanning an array containing coordinate values of points. Succinct data structures are therefore expected to answer queries in sublinear time.

The space complexity of lg *n* + o(lg *n*) bits in the definition of succinct data structures indicates that the size of auxiliary indexing data structures added to the data is negligibly small compared with the size of the data itself (lg *n* bits). In other words, the space complexity of succinct data structures asymptotically matches the information-theoretic lower bound when *n* → ∞.

We compute the information-theoretic lower bound for representing a set of points with integer coordinates. Assume that *i*-th coordinate value takes integer values from 0 to *Ui* − 1. Because the number of grid points is *d*−1 *i*=0 *Ui* , the number of different sets of *n* points is

$$
\begin{pmatrix}
\prod\_{i=0}^{d-1} U\_i \\
\end{pmatrix}.
$$

By using Stirling's approximation formula

$$
\log n! = n \log n - n + \mathcal{O}(\log n),
$$

we obtain

$$\begin{split} \log\binom{U}{n} &= \log U! - \log(U-n)! - \log n! \\ &= U \log U - U - (U-n)\log(U-n) + (U-n) - n\log n + n + \mathcal{O}(\log U) \\ &= U \log\frac{U}{U-n} + n \log \frac{U-n}{n} + \mathcal{O}(\log U) \\ &= U \log\left(1 + \frac{n}{U-n}\right) + n \log \frac{U}{n} \left(1 - \frac{n}{U}\right) + \mathcal{O}(\log U) \\ &= U \left(\frac{n}{U-n} - \Theta\left(\left(\frac{n}{U-n}\right)^2\right)\right) + n \log \frac{U}{n} - \Theta\left(\frac{n^2}{U}\right) + \mathcal{O}(\log U) \\ &= n \log U - n \log n + \Theta\left(n\right). \end{split}$$

Therefore, the information-theoretic lower bound of the size for representing the point set is

$$\lg\left(\prod\_{n=0}^{d-1} \frac{U\_i}{n}\right) = \sum\_{i=0}^{d-1} n \lg U\_i - n \lg n + \Theta\left(n\right).$$

Note that storing coordinate values of the points explicitly using *d*−1 *i*=0 lg *Ui* use *n* lg *n* bit more space than the information-theoretic lower bound.

#### *8.2.2 Assumptions on Point Sets*

Because data structures such as kd-tree or range trees that have linear or larger space complexities usually store the coordinates of points in a plain format, we do not care whether they are integers or real values. However, if we consider succinct data structures, we usually assume that coordinates values are integers from 0 to *n* − 1. We also assume that for any points *p*, *q* ∈ *P* and any *i* ∈ [*d*], the *i*-th coordinate value *pi* of *p* and the *i*-th coordinate value *qi* of *q* are different. Although this assumption may appear to be unrealistic and too strong, for the orthogonal range search problem, it is known that an arbitrary point set on R*d* can be transformed into a point set on [*n*] *d* [7].

Consider a set *P* of *n* points on R*d* . We create another point set *P* on [*n*] *d* as follows. The set *P* also contains *n* points and there is a one-to-one correspondence between points in *P* and points in *P* . Assume that *p* ∈ *P* corresponds to *p* ∈ *P* . Then, the *i*-th coordinate value *p i* of *p* is then defined from the *i*-th coordinate value *pi* of *p* as

#### 8 Orthogonal Range Search Data Structures 127

$$p'\_i = \#\{q \in P \mid q\_i < p\_i\}.\tag{8.1}$$

That is, the *i*-th coordinate value of *p* is the number of points in *P* such that the *i*-th coordinate value is smaller than *pi* . This is called the rank value of *p* with respect to the *i*-th coordinate value, and the transformation is called the transformation into rank space. We use arrays *C*0,*C*1,...,*Cd*−1 each of length *n*. The array *Ci* stores the *i*-th coordinate values of points in *P* in increasing order.

By using the point set *P* on the rank space and the arrays *Ci* (*i* = 0,..., *d* − 1) that contain the original coordinate values of the points in *P*, we can reduce the problem of orthogonal range search on the original point set *P* into that on *P* . Assume that a query range *Q* = - *l* (*Q*) 0 , *u*(*Q*) 0 ×···× - *l* (*Q*) *d*−1, *u*(*Q*) *d*−1 ⊂ R*d* is given for a point set *P*. From the construction of *P* , there exists a range *Q* = - *l* (*Q* ) 0 , *u*(*Q* ) 0 × ···× - *l* (*Q* ) *d*−1 , *u*(*Q* ) *d*−1 ⊂ [*n*] *d* such that

$$p \in \mathcal{Q} \iff p' \in \mathcal{Q'}.$$

The boundaries of this *Q* are computed by

$$\begin{aligned} l\_i^{(\mathcal{Q}')} &= \# \left\{ p \in P \; \middle| \; p\_i < l\_i^{(\mathcal{Q})} \right\} \\ \mu\_i^{(\mathcal{Q}')} &= \# \left\{ p \in P \; \middle| \; p\_i \le \mu\_i^{(\mathcal{Q})} \right\} - 1. \end{aligned}$$

These are computed in O(*d* lg *n*) time by binary searches on the arrays *Ci* . Then, the counting query is performed by using *Q* . For the reporting query, after finding a point *p* ∈ *P* which is included in the query range *Q* in the rank space, we need to recover the original coordinates of the point *p* ∈ *P*. This is done in O(*d*) time using the arrays *Ci* containing the coordinates of the original points by

$$p\_i = C\_i[p'\_i].$$

Thus, an orthogonal range search problem on R*d* can be transformed into that on [*n*] *d* . Note that if coordinates are transformed as in Eq. (8.1), the identical coordinate values in R*d* are transformed into identical coordinate values in [*n*] *d* . By shifting values by one for the identical coordinate values, we can transform the coordinate values so that for any two distinct points *p* , *q* ∈ *P* and any *i* ∈ [*d*], the *i*-th coordinate value *p i* of *p* is different from the *i*-th coordinate value *q i* of *q* .

If the original points have integer coordinate values, we can reduce the space [19]. Consider the case where *P* is a point set on [*U*] *d* , that is, each coordinate value takes an integer value from 0 to *U* − 1. In this case, the point set *P* in the rank space does not change. However, we store the coordinates of the original point set *P* in a different way. We store them using multi-sets *M*0, *M*1,..., *Md*−1, each of which corresponds to one of the *d* dimensions. The multi-set *Mi* stores the *i*-th coordinate value of the points in *P*. We use the data structure of [22] to store multi-sets.

**Lemma 8.1** *There exists a data structure using n* lg(*U*/*n*) + O(*n*) *which supports a selectm query on a multi-set Mi in constant time.*

A selectm query on a multi-set *M* finds the *j*-th smallest element in *M*. That is, *Ci*[*j*] is obtained by finding the *j*-th smallest element in array *Ci* . Therefore, if a query range *Q* on [*U*] *d* is given, it can be transformed into a query range *Q* on the rank space by binary searches using selectm queries, and the original coordinate values are obtained by *d* many selectm queries.

Assume that there exists a succinct data structure *D* for a point set *P* on [*n*] *d* . Then, the space complexity of *D* is(*d* − 1)*n* lg *n* + (*d* − 1) · o(*n* lg *n*) bits, as shown in Sect. 8.2.1. If we add *d* data structures of Lemma 8.1, the total space complexity becomes *dn* lg *U* − *n* lg *n* + (*d* − 1) · o(*n* lg *n*) bits. This is succinct for the point set *P* on [*U*] *d* . Therefore, if there exists a succinct data structure for a point set on [*n*] *d* , we can construct a succinct data structure for a point set on [*U*] *d* . From here onward, we consider only point sets on [*n*] *d* .

#### **8.3 kd-Tree**

kd-tree [1] is a well-known data structure that partitions the space recursively. It is used not only for the orthogonal range search problem, but also for the nearest neighbor search problem.

#### *8.3.1 Construction of kd-Trees*

We explain the algorithm for constructing a kd-tree of a point set *P* for the twodimensional case. First, we find the point *p* for which the *x*-coordinate is the median of the point set *P*, and store *p* at the root of the kd-tree. Next, we divide the set *P* \ {*p*} into two: the set *P*left that stores points with *x*-coordinates smaller than that of *p*, and the set *P*right that stores points with *x* coordinates larger than that of *p*. We add two children vleft, vright to the root of the kd-tree. Next, from *P*left (*P*right), we find *p*left (*p*right) for which the *y*-coordinate is the median of the set, and we store *p*left (*p*right) in vleft (vright). Similarly, we divide the set *P*left \ {*p*left} (*P*right \ {*p*right}) into two subsets according to *y*-coordinates, find medians with respect to *x*-coordinates, and store them in children of vleft (vright), and repeat this recursively. Figure 8.2 shows an example of partitioning a point set.

For a *d*-dimensional space, we partition the space based on the first dimension, the second dimension, and so on. After partitioning the space based on the *d*-th dimension, we use the first dimension again.

**Fig. 8.2** Partitioning of a space based on the point set (left) and the corresponding kd-tree (right). Points *A*, *B*,*C*,... correspond to nodes *a*, *b*, *c*,.... The range corresponding to node *h* is shown in gray in the left figure

#### *8.3.2 Range Search Algorithm*

An important concept for understanding range searches using a kd-tree is the correspondence between nodes of the kd-tree and ranges. In Sect. 8.3.1, we explained that each node of the kd-tree stores a point. We can also consider that each node corresponds to an orthogonal range. Let *V*(v) denote the point in *P* stored in node v and *R*(v) denote the corresponding range. Then *R*(v) is defined as follows:

For example, in Fig. 8.2, the range *R*(*h*) corresponding to node *h* is the gray area.

The algorithm for reporting queries using a kd-tree is as follows. The algorithm searches the space by traversing tree nodes from the root. Each time a node v is visited, the algorithm checks whether the corresponding point *V*(v) (∈ *P*) is contained in the query range *Q* or not. If the range *R*(v) is fully contained in the query range *Q*, the algorithm outputs all the points stored in the sub-tree rooted at v. If *R*(v) and *Q* has no intersection, the algorithm terminates the search of the sub-tree. For a counting query, instead of outputting all the points when *R*(v) is contained in *Q*, the algorithm finds and accumulates the size of the sub-tree rooted at v. Although it may seem impossible to execute the algorithm since the range *R*(v) for node v is not explicitly stored in the kd-tree, if the range *R*(v) for node v is known, then we know the coordinate values of the hyperplane partitioning the range from the coordinate values of point *V*(v), and we can compute *R*(vleft) and *R*(vright). Therefore, we can execute the algorithm by keeping the range *R*(v) during the search.

#### *8.3.3 Complexity Analyses*

The time complexity of kd-trees is analyzed in [13]. A counting query takes O *d* · *n d*−1 *d* + 2*d* time. In general, we assume *d* is a constant and write the complexity as O *n d*−1 *d* . For a reporting query, we output all coordinates of points in *Q*. Because a point can be output in constant time, the query time complexity is O *n d*−1 *d* + *k* .

If *d* ≥ lg *n*, the height of the kd-tree is at most *d*, and therefore the space is partitioned at most *d* times. Then, it is necessary to traverse all the nodes and a query takes O(*n*) time.

#### **8.4 Wavelet Tree**

Wavelet tree is a succinct data structure supporting various queries on strings and integer sequences efficiently. It was originally proposed for representing compressed suffix arrays [9], but it later became known that wavelet tree can support more operations [18]. Orthogonal range search in two-dimensional space is one of these operations [16].

#### *8.4.1 Construction*

The two-dimensional point sets *P* that can be represented directly using wavelet tree are those where the coordinates take integer values from 1 to *n* and the *x*-coordinate values are all distinct. As explained in Sect. 8.2.2, without loss of generality, we can transform any point set into a point set in [*n*] *d* space. For such a two-dimensional point set *P*, consider an integer sequence *C* that contains the *y*-coordinates of the points in increasing order of *x*-coordinates. For example, for the point set in Fig. 8.3, the corresponding integer sequence *C* is 4, 2, 7, 5, 0, 3, 1, 6. For this sequence *C*, we construct a wavelet tree as follows.

First, we consider that the root of the wavelet tree corresponds to *C*. Note that we do not store *C* directly in the wavelet tree. We then focus on the most significant (highest) bit of the lg *n*-bit binary representation of each integer in *C*. If it is 0 (1), the integer is moved into the left (right) child of the root. We consider that each child node of the root corresponds to an integer sequence containing the numbers in the original array *C* in the same order. For example, in the example in Fig. 8.3, integers from 0 to 3 go to the left child, and integers from 4 to 7 go to the right child. Therefore, the left child corresponds to an integer sequence 2, 0, 3, 1, and the right child 4, 7, 5, 6.

**Fig. 8.3** A two-dimensional point set *P* (left) and the corresponding wavelet tree (right)

Next, for each integer sequence of child nodes, we focus on the second most significant bit of the binary representation of each number. We move a number with 0 bit to the left, and a number with 1 bit to the right. Similarly, we repeat this until the integer sequence of a node consists of the identical integer.

Note that we do not store integer sequences in nodes of the wavelet tree. In each node, we store a bit string of the same length as the corresponding integer sequence. The *i*-th bit of the bit string is 0 (1) if the *i*-th integer in the integer sequence goes to the left (right) child. In other words, a bit string stored in a node of depth *l* is the concatenation of the (*l* + 1)-th highest bit of each integer in the integer sequence corresponding to the node. In the example in Fig. 8.3, the integer sequence corresponding to the root node is 4, 2, 7, 5, 0, 3, 1, 6, and because integers from 0 to 3 go to the left child and integers from 4 to 7 go to the right child, the bit string stored in the root node is 1, 0, 1, 1, 0, 0, 0, 1. Note that we do not store bit strings at leaf nodes. We show the information stored in the wavelet tree in the right tree in Fig. 8.3. Only bit strings drawn above the dark gray rectangles, that is, those in the lower row of each node, are stored.

Note that although it may seem impossible to recover the original information (the integer sequence) from these bit strings, it is possible. Consider the recovery of the fourth integer of the wavelet tree in Fig. 8.3 (right). From the bit string stored in the root node, we know that the first bit of the integer is 1. Because this 1 bit corresponding to the fourth integer is the third 1 in the bit string, we know that the integer to be recovered corresponds to the third bit of the bit string in the right child of the root node. If we look at the third bit of the right child, we know that the second bit of the integer is 0. Further, this 0 bit is the second 0 in the bit string, the integer to be recovered corresponds to the second bit of the left child of the current node. Finally, from the second bit of the left child, we know the last bit of the integer to be recovered is 1. Therefore, the fourth integer is 101 in binary, that is, 5. This is shown in Fig. 8.4.

In this recovery operation, we need to compute the number of zeros/ones in the first *i* bits of a bit string. This operation is also used in the range search algorithm in the next section. If we look at bits one by one from the beginning of a bit string, it takes O(*n*) time, which is too slow. We therefore represent the bit string of each node by the following data structure [5, 12].

**Lemma 8.2** *For a bit string of length n, there exists a data structure using n* + o(*n*) *bits which answers a rank/select query in constant time, where the rank query* rank*b* (*B*,*i*) *is to count the number of b bits (b* = 0, 1*) in the bits from B*[0] *to B*[*i*] *(i* ≥ 0*) of a bit string B, and the select query* select*b* (*B*,*i*) *is to return the position of the i -th b (i* ≥ 1*, b* = 0, 1*) in a bit string B.*

The select query is also necessary for range searches using a wavelet tree.

#### *8.4.2 Range Search Algorithm*

We explain how to solve the two-dimensional range search problem using a wavelet tree. First, we explain the counting query, which is performed by a recursive function as in Algorithm 1. For a query range *Q* = [*l*,*r*]×[*b*, *t*], the argument of the function is WTCounting(*l*,*r*, *b*, *t*, vroot, 0, 2lg *n* − 1), where vroot is the root node of the wavelet tree. The left (right) child of node v is represented by vleft (vright). The bit string stored in node v is represented by v.*B*.

We explain the algorithm in Fig. 8.5 using the example of searching a range *Q* = [1, 6]×[1, 4] for the point set *P* in Fig. 8.3.

The search algorithm traverses the tree from the root. During the search, the algorithm keeps the interval *I* of an integer sequence (or bit string) corresponding to an interval of the *x*-coordinate of the query range. In the example in Fig. 8.5, we focus on the interval *I* = [1, 6] at the root node. To move to the left child, we need to compute the interval corresponding to the query range. This is done by a rank query that counts the number of zeros from the beginning of the bit string to a specified position. In the bit string stored in the root node, the number of zeros from the beginning to the 0-th position (in general, if the interval is *I* = [*l*,*r*], to (*l* − 1)-th

**Fig. 8.5** Behavior of the algorithm when searching a range of[1, 6]×[1, 4]for the two-dimensional point set in Fig. 8.3

**Algorithm 1** WTCounting(*x*1, *x*2, *y*1, *y*2, v, *a*, *b*)

**Input:** A node v of the wavelet tree and an interval [*x*1, *x*2] in the corresponding bit string, the interval [*a*, *b*] of *y* coordinate corresponding to node v, and the interval [*y*1, *y*2] of *y* coordinate for the query range.

**Output:** The number of points stored in the sub-tree rooted at v and contained in *Q*.

```
1: if x1 > x2 then
2: return 0
3: else if [a, b]∩[y1, y2] = ∅ then
4: return 0
5: else if [a, b]⊆[y1, y2] then
6: return x2 − x1 + 1
7: end if
8: xl
    1 ← rank0 (v.B, x1 − 1)
9: xl
    2 ← rank0 (v.B, x2) − 1
10: xr
     1 ← x1 − xl
                 1
11: xr
     2 ← x2 − xl
                 2 − 1
12: m ← (a + b)/2
13: return WTCounting(xl
                            1, xl
                                2, y1, y2, vleft, a, m)
           +WTCounting(xr
                             1 , xr
                                 2 , y1, y2, vright, m + 1, b)
```
position) is 0, so we know the interval corresponding to the query starts at position 0. Because the number of zeros from the beginning to the 6-th position (in general, if the interval is *I* = [*l*,*r*], to *r*-th position) is four, we know the interval ends at position 3. Thus, we obtain the interval *I* = [0, 3] for the left child. Similarly, for the right child, by using rank queries counting the number of ones, we can obtain the interval *I* = [1, 2].

We repeat this process by going down the tree maintaining an interval. When we reach a leaf, we can determine if the *y*-coordinate of the point is included in the query range. However, we can sometimes determine this at an earlier stage. For example, in Fig. 8.5, after obtaining the interval *I* = [1, 2] at the left child of the root, for the right child of the current node the interval of the *y*-coordinate corresponding to the node is [2, 3], which is completely included in the interval [1, 4] of the *y*-coordinate of the query range. Therefore, for the two points we focus on at this node, both the *x*- and *y*-coordinates are included in the query range, and we found two points in the query range. However, after computing the interval *I* = [1, 2] for the right child of the root, the interval of the *y*-coordinate corresponding to the right child of the current node is [6, 7], which has no intersection with the interval [1, 4] of the *y*-coordinate of the query range. We do not need to further search the sub-tree.

As observed above, in a range search using a wavelet tree, if the query range is *Q* = [*l*,*r*]×[*b*, *t*], we first focus on points for which the *x*-coordinates are contained in *Q*, that is, contained in the range [*l*,*r*]×[0, *n* − 1]. Next, the process of traversing down the tree corresponds to partitioning the range into two according to the *y*coordinate. If an obtained range is completely contained in the query range, or does not intersect with the query range, we terminate searching the sub-tree.

For counting queries, it is sufficient to sum the number of points. For reporting queries, the extra work of computing the coordinates of the points is also required. This is shown in Algorithm 2.

The outline of the reporting query is the same as the counting query. In Algorithm 1, we obtain the number of points in Line 2. We change it one by one to output coordinates of points corresponding to the interval [*x*1, *x*2] of the bit string v.*B*. The *x*- and *y*-coordinates of each point are obtained by WTReportX

**Input:** A node v of the wavelet tree, the interval [*x*1, *x*2] of the bit string stored in it, the interval [*a*, *b*] of *y* coordinates corresponding to the range for v, and the interval [*y*1, *y*2] of *y* coordinates for the query range. **Output:** Coordinates of point stored in the sub-tree rooted at v and contained in *Q*. 1: **if** *x*1 > *x*2 **then** 2: terminate 3: **else if** [*a*, *b*]∩[*y*1, *y*2] = ∅ **then** 4: terminate 5: **else if** [*a*, *b*]⊆[*y*1, *y*2] **then** 6: **for** *i* = *x*1 to *x*2 **do** 7: *x* ← WTReportX(v,*i*) 8: *y* ← WTReportY(v,*i*, *a*, *b*) 9: Output (*x*, *y*) 10: **end for** 11: **end if** 12: *xl* 1 ← rank0 (v.*B*, *x*1 − 1) 13: *xl* 2 ← rank0 (v.*B*, *x*2) − 1 14: *xr* 1 ← *x*1 − *xl* 1 15: *xr* 2 ← *x*2 − *xl* 2 − 1 16: *m* ← (*a* + *b*)/2 17: WTReporting(*xl* 1, *xl* 2, *y*1, *y*2, vleft, *a*, *m*) 18: WTReporting(*xr* 1 , *xr*

2 , *y*1, *y*2, vright, *m* + 1, *b*)

# **Algorithm 3** WTReportX(v,*i*)

**Input:** A node v of the wavelet tree and an integer *i*.

**Output:** The *x* coordinate value of the point corresponding to the *i*-th bit of the bit string stored in v.

```
1: if v is the root then
2: return i
3: else if v is the left child of vparent then
4: i ← select0

                  vparent.B,i + 1

5: return WTReportX(vparent,i)
6: else
7: i ← select1

                  vparent.B,i + 1

8: return WTReportX(vparent,i)
9: end if
```
# **Algorithm 4** WTReportY(v,*i*, *a*, *b*)

**Input:** A node v of the wavelet tree, the interval [*a*, *b*] of *y* coordinate corresponding to the range for v, and an integer *i*.

**Output:** The *y* coordinate value of the point corresponding to the *i*-th bit of the bit string stored in v.

```
1: if a = b then
2: return a
3: else if v.B[i] = 0 then
4: i ← rank0 (v.B,i) − 1
5: return WTReportY(vleft,i, a, (a + b)/2)
6: else
7: i ← rank1 (v.B,i) − 1
8: return WTReportY(vright,i, (a + b)/2 + 1, b)
9: end if
```
and WTReportY, respectively. The algorithm WTReportY for computing the *y*coordinate (Algorithm 4) is similar to the algorithm for recovering a value of the original integer array explained in Sect. 8.4.1. We compute the *y*-coordinates by traversing down the tree using rank queries.

In contrast, the algorithm WTReportX for computing the *x*-coordinate (Algorithm 3) traverses up the tree using select queries. We explain this by example. In Fig. 8.5, assume that at node v, which is the right child of the left child of the root, we find that points corresponding to the interval *I* = [0, 1] are contained in the query range. Consider the computation of the *x*-coordinate of the point corresponding to the bit v.*B*[1]. First, the node v we focus on is the right child of its parent. We find the position of the second 1 in the parent by a select query. Then we know that the point corresponds to the bit v .*B*[2] in the parent node v . Next, because the current node is the left child of the parent (the root), we find the position of the third 0 in the bit string of the parent by a select query. Now we know that the point corresponds to the bit *r*.*B*[5] at the root node *r*. That is, the *x*-coordinate of the point is 5.

As shown above, we can traverse the nodes of the wavelet tree using rank and select queries on bit strings. For range searches, we traverse down the tree from the root computing the intervals of the *x*-coordinate corresponding to the query range. If we find a node where the corresponding interval of the *y*-coordinate is contained in the query range, we answer the query by computing the length of the interval or coordinate values by traversing the tree.

#### *8.4.3 Complexity Analyses*

We now analyze the space complexity of the wavelet tree and query time complexities for the orthogonal range search problem.

First, we analyze the space complexity. The height of the wavelet tree is lg *n*. The total length of bit strings stored in the nodes with the same depth is always *n*. Therefore, the total length of all the bit strings in the wavelet tree is *n* lg *n*. We can concatenate all the bit strings and store only a long bit string. Then it is not necessary to store the tree structure of the wavelet tree. By using the data structure of Lemma 8.2 for this long bit string, the space complexity is *n* lg *n* + o(*n* lg *n*) bits in total.

Next, we consider query time complexities. For a counting query, we consider the number of visited nodes. In the wavelet tree, each time we traverse an edge toward a leaf, points with small *y*-coordinates go to the left child, and points with large *y*coordinates go to the right child. At leaves we can consider that all points are sorted in increasing order of *y*-coordinates. This means that leaf nodes corresponding to the interval of *y*-coordinates of the query range exist in a consecutive place in the wavelet tree. Now, consider the set *M* of nodes of the wavelet tree defined as follows. The set *M* contains a maximal node v such that the *y*-coordinates corresponding to the leaf nodes in the sub-tree rooted at v are contained in the query range, that is, the *y*-coordinates of the leaves in the sub-tree of v are contained in the query range but the sub-tree of the parent of v contains some node for which the corresponding *y*-coordinate is not contained in the query range. This is the set of nodes from which we do not further search the sub-tree for a counting query using the wavelet tree, and in Fig. 8.6, it is shown as dark gray nodes.

Let *A* be the set of nodes that are ancestors of nodes of *M*. This is the set of nodes visited before reaching nodes of *M* which are shown as light gray nodes in Fig. 8.6. The number of nodes visited in a counting query is then |*A*|+|*M*|. We now consider the size of *M* and *A*.

For the size of the set *M*, the following lemma holds.

## **Lemma 8.3** *It holds* |*M*| = O(lg *n*)*.*

*Proof* (Lemma 8.3) The set *M* is constructed as follows. Let *M* be the set of leaf nodes of the wavelet tree corresponding to the interval of *y*-coordinates in the query range. For the nodes of *M* , if two nodes v1 and v2 have a common parent node v, we remove v1 and v2 from *M* and add v to *M* . By repeating this process until there are no such pairs of nodes, the set *M* coincides with *M*.

For each depth of the wavelet tree, the number of nodes of depth belonging to *M* is then at most two, because if there exist more than two nodes, two of them must have the same parent. This completes the proof that |*M*| = O(lg *n*).

**Fig. 8.6** Nodes visited by a counting query. We traverse light gray nodes, and when we reach a dark gray node, we do not further search the nodes below it

For the size of the set *A*, the following lemma holds.

**Lemma 8.4** *It holds* |*A*| = O(lg *n*)*.*

*Proof* (Lemma 8.4) Consider a node v in the set *A*. In the set of leaf nodes in the subtree rooted at v, there must exist a leaf node where the corresponding *y*-coordinate is included in the query range and a leaf node where the corresponding *y*-coordinate is not included in the query range. Therefore, for each depth of the wavelet tree, there are at most two such nodes in *A*, because if there exists more than two such nodes, for a node in the middle, the corresponding *y*-coordinates of the leaves in the sub-tree rooted at that node are contained in the query range. This completes the proof that |*A*| = O(lg *n*).

From the above discussion, the number of nodes visited in a counting query is |*A*|+|*M*| = O(lg *n*). When we visit a new node, we use a constant number of rank queries. Because a rank query takes constant time (Lemma 8.2), the time complexity of a counting query using the wavelet tree is O(lg *n*).

For a reporting query, it is necessary to compute coordinates of points in the query range. As explained in Sect. 8.4.2, *x*-coordinates are computed by traversing up the tree and *y*-coordinates are computed by traversing down the tree, with the coordinates of each point computed by visiting O(lg *n*) nodes. Moving to an adjacent node in the wavelet tree is done by a constant number of rank/select queries, and each rank/select query takes constant time (Lemma 8.2). Therefore, the coordinates of a point are obtained in O(lg *n*) time, and the time complexity for a reporting query using the wavelet tree is O((1 + *k*)lg *n*), where *k* is the number of output points.

We obtain the following theorem.

**Theorem 8.1** *The space complexity of the wavelet tree representing a twodimensional point set on* [*n*] 2 *is n* lg *n* + o(*n* lg *n*) *bits, and a counting query takes* O(lg *n*)*time, and a reporting query takes* O((*k* + 1)lg *n*)*time, where k is the number of points to enumerate.*

As shown in Sect. 8.2.1, the information-theoretic lower bound for a point set on [*n*] 2 is *n* lg *n* + O(*n*) bits. Therefore, the wavelet tree is a succinct data structure.

**Theorem 8.2** *Let P be a set of points on M* = [1..*n*]×[1..*n*]*in which all points have distinct x -coordinates. Then, there exists a data structure using n* lg *n* + o(*n* lg *n*) *bits that answers a counting query in* O(lg lg *n*) *time and a reporting query in* O((1 + *k*)lg *n*/ lg lg *n*) *time, where k is the number of points to output.*

## **8.5 Proposed Data Structure 1: Improved Query Time Complexity**

This data structure uses the idea of adding data structures to the kd-tree to improve the query time complexity [20]. First, we explain the idea of [20] in Sect. 8.5.1. Next, we explain the algorithm of range search in Sect. 8.5.3, and analyze the time complexity in Sect. 8.5.4.

#### *8.5.1 Idea for Improving the Time Complexity of the kd-Tree*

The method proposed in [20] improves the query time complexity of the kd-tree by adding *d* many wavelet trees to the kd-tree such that the term *n*(*d*−1)/*d* is replaced by *n*(*d*−2)/*d* (lg *n* if *d* = 2), at the cost of increasing the total complexity by a factor of O(lg *n*). Note that we assume point sets are on [*n*] *d* .

First, we construct the kd-tree for a given set *P* of points in the *d*-dimensional space. Next, we label the nodes of the kd-tree with numbers based on the inorder traversal of a binary tree defined as follows:

Figure 8.7 shows an example of a point set (left) and numbers assigned based on the inorder traversal of the kd-tree of the set (right).

Next, we make point sets *Pi* (*i* = 0,..., *d* − 1) with *n* points on [*n*] 2. The two-dimensional point set *Pi* is created as follows. If a point *p* in the original *d*dimensional point set *P* has the *i*-th coordinate value *pi* and the inorder position of the node of the kd-tree containing *p* is *j*, we add point (*j*, *pi*) to *Pi* . Figure 8.8 shows the point sets *P*0, *P*1 created from the point set in Fig. 8.7.

**Fig. 8.7** A two-dimensional point set (left) and the corresponding kd-tree (right). The numbers of nodes are assigned by an inorder traversal of the kd-tree. The dashed lines in the left figure show the partition of the space by the kd-tree

**Fig. 8.8** Two-dimensional point sets obtained from the point set in Fig. 8.7

From these two-dimensional point sets *P*0,..., *Pd*−1, we construct wavelet trees *W*0,..., *Wd*−1. The wavelet trees *Wi* can be thought of as constructed from an integer sequence *Ai* containing the *i*-th coordinate value of points in *P* in the order of the kd-tree.

These data structures can be used for range searches as follows. Given a query range *Q*, we perform the original search using the kd-tree. In the original algorithm, as explained in Sect. 8.3, we traverse the kd-tree and shrink the range *R*(v), and when CDeg(*R*(v), *Q*) = *d* (i.e., *R*(v) ⊆ *Q*), we know that all the points in the sub-tree rooted at v are contained in *Q*. By using the *d* wavelet tree, we can terminate the search when CDeg(*R*(v), *Q*) = *d* − 1. Assume that when a node v is visited, *R*(v) is contained in *Q* for all dimensions except for *i*. The inorder numbers of nodes in the sub-tree rooted at node v have consecutive values. Let [*a*, *b*] be the interval for the numbers. Then, the points in this interval are contained in *Q* except for dim. *i*. This implies that points in *Pi* that are contained in the range [*a*, *b*]×[*l* (*Q*) *i* , *u*(*Q*) *i* ] are contained in *Q* even for dim. *i*. Therefore, after finding the node v, it is sufficient to search the range [*a*, *b*]×[*l* (*Q*) *i* , *u*(*Q*) *i* ] of *Pi* using wavelet trees *Wi* .

The number of nodes of the kd-tree visited by this method is O *n*(*d*−2)/*d* (O(lg *n*) for the case *d* = 2). The search of the last dimension using the wavelet tree takes O(lg *n*) time for a counting query. Therefore, the time complexity for a counting query using the kd-tree is improved to O *n*(*d*−2)/*d* lg *n* (O lg2 *n* for the case *d* = 2).

#### *8.5.2 Index Construction*

We now explain the proposed data structure. First, we construct the kd-tree for a given point set *P*. Note that this kd-tree is temporarily built in order to construct our data structure, and is not included in the final structure. Next, as in Sect. 8.5.1, we number the nodes of the kd-tree by an inorder traversal, and create *d* many two-dimensional point sets *P*0,..., *Pd*−1. For each *Pi* , we create the data structure of [3]. Let *Bi* be this data structure. Finally, we discard the kd-tree. The final data structure consists of *B*0,..., *Bd*−1.

#### *8.5.3 Range Search Algorithm*

We explain the algorithm for a reporting query using the data structure explained in the previous section. The pseudocode is shown in Algorithm 5. This algorithm simulates a search of the kd-tree using *B*0,..., *Bd*−1. We explain it in comparison with the search algorithm of the kd-tree. Note that we explain the algorithm assuming the inorder number of each node v of the kd-tree is also assigned to the point *V*(v) stored in v. That is, if we say a point with number *j*, it is the point stored in the node with inorder position *j*. We also assume that for an interval [*a*, *b*] of point numbers, *R*([*a*, *b*]) denotes the range containing points that have numbers in [*a*, *b*]. In Algorithm 5, the interval[*a*, *b*] of point numbers always corresponds to the interval of inorder numbers of nodes in the sub-tree rooted at a node v of the kd-tree. Therefore, *R*([*a*, *b*]) coincides with *R*(v).

If we use the kd-tree, we shrink the focused range *R*(v) by going down the tree. In the proposed method, by shrinking the interval [*a*, *b*] of point numbers, we reduce the corresponding range *R*([*a*, *b*]). Because the kd-tree stores the point *V*(v) corresponding to a node v, we can obtain the information of the point used for

# **Algorithm 5** Reporting([*a*, *b*], *Q*)

**Input:** An interval of point numbers [*a*, *b*] and a query range *Q*.

**Output:** Points with numbers in [*a*, *b*] and which are contained in range *Q*.

1: **if** *Deg*(*R*([*a*, *b*]), *Q*) = *d* − 1 **then** 2: For the last dimension - *i* such that *R*([*a*, *b*]) is not yet contained in *Q*, search [*a*, *b*] × *l* (*Q*) *i* , *u*(*Q*) *i* of *Pi* using *Bi* and enumerate points contained in *Q*. For each point, compute coordinates using *B*0,..., *Bd*−1. 3: **else if** *R*([*a*, *b*]) has no intersection with *Q* **then** 4: terminate 5: **else if** *a* = *b* **then** 6: Examine the point with number *a*. If it is in *Q*, output it. 7: **else** 8: *c* ← (*a* + *b*)/2 9: Output the point with number *c* if it is in *Q*. 10: Reporting([*a*, *c* − 1], *Q*) 11: Reporting([*c* + 1, *b*], *Q*) 12: **end if**

partitioning the space. In contrast, in the proposed method, points are not explicitly stored. However, if the focused interval [*a*, *b*] coincides with the interval of inorder numbers for the sub-tree rooted at a node v, we find *c* = (*a* + *b*)/2 is the number of the points used for partitioning.1 Furthermore, the intervals [*a*, *c* − 1] and [*c* + 1, *b*] correspond to the intervals of the numbers for sub-trees rooted at the left and right child of v, respectively. Therefore, by a recursive search of Algorithm 5, we can obtain the correct partitioning points.

For the range *R*([*a*, *b*]), we can compute the ranges after a partition from the range before partition and the coordinates of the point used for partitioning similarly to the case of the kd-tree.

#### *8.5.4 Complexity Analyses*

We now analyze the complexities of the algorithm. First, we consider its space complexity.We use *d* data structures of Bose et al. [3] each of which uses *n* lg *n* + o(*n* lg *n*) bits as in Theorem 8.2. The total space complexity is then *dn* lg *n* + o(*dn* lg *n*) bits.

Next, we consider the query time complexity. If we use the same analysis as in [20], assuming *d* is a constant, we can show the number of nodes corresponding to cells with containment degree of at most *d* − 1 is O *n d*−2 *d* . Here, we derive the query time complexity using a novel analysis for non-constant *d*.

1 In the kd-tree, at each depth, we partition the space by the median of the point set with respect to a dimension, and therefore *c* = (*a* + *b*)/2 is the number of the point used for partitioning. If the point set contains an even number of points, we can obtain the correct partitioning point using a predetermined rule.

The proposed method partitions the space for each dimension in order, in the same fashion as for the kd-tree. As in [20], we define a series of partitions with respect to dim. 0 to dim. d-1 as a cycle. We then calculate the number *Tm* (*n*, *d*) of nodes at which the containment degree with respect to *Q* is at most *d* − 2 in the *m*-th cycle. When the (*m* − 1)-th cycle has finished, the space is partitioned into 2*d*(*m*−1) many cells. Among them, we count the number of cells for which the containment degree with respect to *Q* is at most *d* − 2. These cells contain a (*d* − 2)-dimensional face of *Q* (an edge of a cuboid if *d* = 3). A (*d* − 2)-dimensional face of a *d*-dimensional orthogonal range *Q* is obtained by choosing two dimensions from the *d*-dimensions and choosing the upper side or the lower side of the range for each of the two dimensions. Therefore, *Q* has *d* 2 22 many (*d* − 2)-dimensional faces. When the (*m* − 1)-th cycle has finished, because each dimension is partitioned into 2*m*−1 cells, the number of cells containing a (*d* − 2)-dimensional face is at most 2(*m*−1)(*d*−2) . Then after the (*m* − 1)-th cycle, the number of cells to be searched is at most

$$
\binom{d}{2} 2^2 \cdot 2^{(m-1)(d-2)}\dots
$$

In the sub-trees rooted at these nodes, the number of nodes in the *m*-th cycle is 2*d* − 1. Therefore, it holds that

$$\begin{aligned} T\_m(n,d) &\le \left(2^d - 1\right) \binom{d}{2} 2^2 \cdot 2^{(m-1)(d-2)} \\ &< 2^3 d (d-1) 2^{(d-2)m} .\end{aligned}$$

Let *N*(*n*, *d*) be the number of nodes for which the containment degree with respect to *Q* is at most *d* − 2. It then holds that

$$\begin{aligned} N(n,d) &= \sum\_{m=1}^{\frac{1}{d}\lg n} T\_m(n,d) \\ &< 2^3 d(d-1) \sum\_{m=1}^{\frac{1}{d}\lg n} 2^{(d-2)m} \\ &= 2^3 d(d-1) \frac{2^{d-2} \left(2^{\frac{d-2}{d}\lg n} - 1\right)}{2^{d-2} - 1} \\ &= \mathcal{O}\left(d^2 \cdot n^{\frac{d-2}{d}}\right). \end{aligned}$$

We use the fact that the containment degree is weakly increasing as we traverse down the tree. In the proposed method, we terminate the search when the containment degree reaches *d* − 1. The visited nodes are then those with containment degree of at most *d* − 2 and their child nodes. There are at most 2*N*(*n*, *d*) such nodes. The proposed method virtually traverses the kd-tree. It takes O *d* lg *n* lg lg *n* time to compute the coordinates of a point stored in a node. When a node for which the containment degree with respect to *Q* is *d* − 1, we search the last one dimension in O lg *n* lg lg *n* time. The time complexity of a counting query is then O *d*3*n d*−2 *d* lg *n* lg lg *n* . For a reporting query, it takes O *d* lg *n* lg lg *n* time to compute the coordinates of a point. The total time complexity is then O*d*3*n d*−2 *d* + *dk* lg *n* lg lg *n* , where *k* is the number of points in *Q*. In summary, we obtain the following:

**Theorem 8.3** *For an orthogonal range search problem on the* [*n*] *d space, there exists a data structure that has space complexity of dn* lg *n* + o(*dn* lg *n*) *bits and which answers a counting query in* O *d*3*n d*−2 *d* lg *n* lg lg *n time and a reporting query in* O *d*3*n d*−2 *d* + *dk* lg *n* lg lg *n , where k is the number of points in the query range.*

## **8.6 Proposed Data Structure 2: Succinct and Practically Fast**

The second proposed method is a data structure that is succinct and practically fast. In this method, we use *d* − 1 many wavelet trees to represent a point set on [*n*] *d* . In Sect. 8.6.1, we explain how to construct the data structure. In Sect. 8.6.2, we explain the algorithm for the orthogonal range search problem. In Sect. 8.6.3, we analyze the space and time complexities.

#### *8.6.1 Index Construction*

In this method, we assume that the points of *P* have distinct values in the 0-th coordinate value.

First, we create length-*n* integer arrays *A*1,..., *Ad*−1. The array *Ai* corresponds to dim. *i*, and stores the *i*-th coordinate value of the points in increasing order of the 0-th coordinate value. Next, for those arrays we create wavelet trees *W*1,..., *Wd*−1. The wavelet trees *Wi* can be considered to represent the two-dimensional point set *Pi* generated from the *d*-dimensional point set *P* by projecting the points onto the plane spanned by the 0-th axis and the *i*-th axis. Figure 8.9 shows an example.

# **Algorithm 6** Report(*Q*)

**Input:** A query range *Q* = - *l* (*Q*) 0 , *u*(*Q*) 0 ×···× - *l* (*Q*) *d*−1, *u*(*Q*) *d*−1 . **Output:** The coordinates of points of *P* contained in *Q*. 1: *D* := ∅ 2: **for** *i* = 1 to *d* − 1 **do** 3: **if** - *l* (*Q*) *i* , *u*(*Q*) *i* [0, *n* − 1] **then** 4: *D* = *D* ∪ {*i*} 5: *ci* := Count *Pi*, - *l* (*Q*) 0 , *u*(*Q*) 0 × - *l* (*Q*) *i* , *u*(*Q*) *i* 6: **end if** 7: **end for** 8: Sort elements *i*1,...,*i*|*D*| of *D* in increasing order of *ci* . 9: *A* := ReportX *Pi*1 , - *l* (*Q*) 0 , *u*(*Q*) 0 × - *l* (*Q*) *i*1 , *u*(*Q*) *i*1 10: **for** *i* = *i*2 to *i*|*D*| **do** 11: **for all** *a* ∈ *A* **do** 12: **if** The *i*-th coordinate of a point for which the 0-th coordinate is *a* and is not contained in - *l* (*Q*) *i* , *u*(*Q*) *i* **then** 13: *A* = *A* \ {*a*} 14: **end if** 15: **end for** 16: **end for** 17: **for all** *a* ∈ *A* **do** 18: Obtain the coordinates of a point for which the 0-th coordinate is *a* and output them. 19: **end for**

#### *8.6.2 Range Search Algorithm*

Next, we explain how to solve the orthogonal range search problem using the data structure (Algorithm 6).

Assume that a query range *Q* = - *l* (*Q*) 0 , *u*(*Q*) 0 ×···× - *l* (*Q*) *d*−1, *u*(*Q*) *d*−1 is given. For each *i* = 1,..., *d* − 1 such that - *l* (*Q*) *i* , *u*(*Q*) *i* = [0, *n* − 1], that is, the dimension *i* used for the search, we count the number of points of *Pi* that are contained in range - *l* (*Q*) 0 , *u*(*Q*) 0 × - *l* (*Q*) *i* , *u*(*Q*) *i* using wavelet trees *Wi* (counting query). Let *m* (= |*D*|) be the number of *i* (= 1,..., *d* − 1) such that - *l* (*Q*) *i* , *u*(*Q*) *i* = [0, *n* − 1], and let *i*1,...,*im* be the sorted ones in increasing order of the number of answers of counting queries.

Using wavelet trees *Wi*1 , we then enumerate only the *x*-coordinates of points of *Pi*1 contained in - *l* (*Q*) 0 , *u*(*Q*) 0 × - *l* (*Q*) *i*1 , *u*(*Q*) *i*1 and store them in a set *A*. For each element *a* of *A* and for each *i* = *i*2,...,*im*, we check whether the *i*-th coordinate of a point for which the 0-th coordinate is *a* is contained in the query range. The elements remaining in *A* correspond to points in the query range. The answer to a counting query is the cardinality of *A*. For a reporting query, we compute coordinates of the points and output them.

The reason we compute the number of points contained in each dimension by a counting query is twofold. Firstly, the *x*-coordinate (the 0-th coordinate) of points contained in the query range with respect to the *i*1-th (and the 0-th) dimension can be output quickly at line 9 of the algorithm if the number of points to enumerate is small. Secondly, in the double loops from line 10 to line 16, we want to reduce the size of *A* as soon as possible.

#### *8.6.3 Complexity Analyses*

Consider the space and time complexities of the proposed method.

For the space complexity, we use *d* − 1 many wavelet trees. Therefore, the space complexity is (*d* − 1)lg *n* + (*d* − 1) · o(lg *n*) bits.

For the query time complexity, let *m* be the number of wavelet trees used in a search. The time to perform *m* counting queries on wavelet tree is O(*m* lg *n*). We then sort *m* integers in O(*m* lg *m*) time. Next, we enumerate the *x*-coordinates of points contained in the query range for the dimension with the minimum number of points. Let *ci*1 = *c*min be the number of points to enumerate. This takes O((1 + *c*min)lg *n*) time. The time to check whether these points are contained in the query range for other dimensions is O((*m* − 1)*c*min lg *n*). Let *d* be the number of dimensions used in the query, then it holds that *m* ≤ *d* . Therefore, the query time complexity can be written as O *d c*min lg *n* + *d* lg *d* .

#### **8.7 Conclusion**

In this chapter, we first reviewed data structures for high-dimensional orthogonal range search. We then proposed two data structures for the problem.

The first one simulates the search of the kd-tree using *d* succinct data structures for two-dimensional orthogonal range search data structures [3]. We improved the query time complexity of KDW-tree while keeping the same space complexity.

The second one is succinct and practically fast. The space complexity is (*d* − 1)*n* lg *n* + (*d* − 1) · o(*n* lg *n*), which is succinct. The worst-case query time complexity is O(*dn* lg *n*), which is not good. However, if the number *d* of dimensions is large but the number *d* of dimensions used in a search is small, it runs fast in practice.

#### **References**

**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

# **Chapter 9 Enhanced RAM Simulation in Succinct Space**

**Taku Onodera**

**Abstract** We describe two recent results on space-efficient functional random access memory (RAM), which is RAM with non-standard functionalities. The first is about oblivious RAM, which enables a remote database to be accessed without revealing to the database owner which part of the database is being accessed. The other is about wear leveling, which enables the number of updates to be balanced among all the memory cells regardless of the content of the computation being performed on the memory.

## **9.1 Introduction**

Random access memory (RAM) underlies most modern computers, and improvements to the RAM itself can have a positive impact on a wide range of applications. For example, faster RAM access makes all RAM-based computations correspondingly faster. Some types of RAM improvements are not just about efficiency but also about *functionality*. An example is virtual memory in operating systems, which enables, among other things, application programs to utilize the memory without concern about cumbersome management issues such as allocation. Generally speaking, this type of RAM improvement functions by using conventional RAM to simulate "enhanced" RAM while introducing some performance overhead.

In this chapter, we describe two such enhanced RAM simulations—*oblivious RAM* (ORAM) and *wear leveling*—with the emphasis on how to minimize the space overhead. These topics were chosen mainly because the authors' knowledge of them, although there are also some conceptual similarities between ORAM and wear

T. Onodera (B)

The University of Tokyo, Tokyo, Japan e-mail<EMAIL_ADDRESS>
© The Author(s) 2022 N. Katoh et al. (eds.), *Sublinear Computation Paradigm*, https://doi.org/10.1007/978-981-16-4095-7\_9

leveling. Other functionality-enhanced RAM simulations include initializable array [2], memory checking [3], locally decodable code [11], and huge random object [8].1

#### **9.2 Oblivious RAM**

#### *9.2.1 Problem*

Suppose you want to outsource a database, stored in RAM, to a server and want to access it in a privacy-preserving way. Although you can hide the data content by encryption, the server can still see which part of the RAM you are accessing. This is a serious issue in the current era of cloud computing. The same problem also appears when one wants to hide the details of software implemented in a physically secure processor that accesses insecure main memory.

Oblivious RAM (ORAM) is the formalization and corresponding solution of this problem. Typically, it works by storing the RAM into some data structure on the server and moves the RAM cells dynamically in the data structure as the user accesses the RAM.

As an example, consider a scheme where the server stores the RAM as-is except that each cell is encrypted by the user's key. To access the *i*th cell, the user performs the following procedure for *j* = 1 to *N* where *N* is the number of cells:

	- For read access, copy the decrypted value to local memory.
	- For write access, change the decrypted value to the new value.

We assume semantically secure encryption when encryption is used in this chapter. In particular, there is an overwhelmingly high probability that the re-encrypted ciphertext looks totally different to the server from the ciphertext before re-encryption regardless of whether the plaintext is updated or not. Thus, no matter what actual access is performed, all that the server can see is that random-looking encrypted cells are updated to still random-looking re-encrypted cells in a fixed scan order. Of course, the access overhead of this method is very large since each cell access takes

1 A huge random object, in this context, is a succinct representation of a pseudorandom object that supports certain queries. For example, a pseudorandom function can be thought of as a huge pseudorandom bitstring that is implicitly represented by a tiny seed and supports efficient random access. This is not a data structure in the conventional sense because the represented object is a pseudorandom bitstring instead of "data."

time linear to the entire RAM size. The purpose of this example is merely to illustrate the kind of security we want to achieve.

We now give a more formal problem description. We have three parties: the user, the server, and the simulator. The simulator models a program that runs in the local environment of the user. The simulator provides the user with an access interface to RAM that we call the virtual RAM while the server provides the simulator with an access interface to RAM that we call the physical RAM. That is, the user gives the simulator a series of queries of the form (type,*i*, *v*) where type ∈ {read, write}, *i* ∈ [*N*], and *v* ∈ {0, 1}*B*. We call these virtual queries. The parameter *N* specifies the number of virtual cells—cells in the virtual RAM—while *B* specifies the size of each virtual cell. Given a virtual query, the simulator gives the server another series of queries of the form (type,*i*, *v*), where type ∈ {read, write}, *i* ∈ [*N* ], and *v* ∈ {0, 1}*B* . We call these physical queries. The simulator, and thus the physical queries, is probabilistic in general.2 The server responds to physical queries in the obvious way. That is, for (read,*i*, ∗) where "∗" means that the third component is arbitrary, the server returns the value of the *i*th physical cell, and for (write,*i*, *v*), the server updates the value of the *i*th physical cell to *v*. If the virtual query from the user is of the form (read,*i*, ∗), the simulator derives the value of the *i*th virtual cell through the interaction with the server and returns it to the user. If the virtual query is of the form (write,*i*, *v*), the simulator updates the value of the *i*th virtual cell to *v*. The simulator must respond to the virtual queries online. We call the sequence of second components of the virtual queries (resp. physical queries) a virtual access pattern (resp. physical access pattern). For a virtual query sequence *q*, let *a*(*q*) denote the physical access pattern induced by *q*. Recall that *a*(*q*) is a random variable in general. The ORAM scheme is secure if *a*(*q*1) and *a*(*q*2) are indistinguishable for any virtual query sequences *q*1 and *q*2 of the same length. There are some variations in the exact meaning "indistinguishable". Typically used meanings of indistinguishability in descending order of security are a) equally distributed, b) statistically closely distributed, and c) computationally indistinguishable. The main performance metric of ORAM includes access overhead, which is the number of physical queries processed for each virtual query, the simulator local space size, and the server space size, which is *B N* bits.

As mentioned above, the simulator models a program running in the local environment of the user. Thus, in practice, we do not distinguish the user and the simulator. For example, we refer to the simulator local space as user space.

The ORAM problem is non-trivial only if the user space is smaller than *BN* bits since otherwise, the simulator can store the entire RAM locally and ignore the server.

2 The simulator of the scan-based method is deterministic, and it is not hard to see that its linear access overhead is optimal if we restrict the simulator to be deterministic. Thus, all simulators of interest are indeed probabilistic.

**Table 9.1** Summary of existing results. † means amortized bound. ≈ log *N* means *O*( *f* (*N*)) for any *f* ∈ ω(log *N*). The constant factor of the user space of [26] is 1

#### *9.2.2 Existing Results*

Table 9.1 gives a summary of some of the existing results. Every method has physical cell size *B* = *B* + (log *N*). There is an (log *N*) lower bound for the access overhead if the user space is at most *N*1− for constant > 0 [14].

There are mainly two types of techniques that are actively studied: hierarchical approaches and tree-based approaches.3 Asymptotically, the state-of-the-art hierarchical method [12] has access overhead matching the lower bound mentioned above while the state-of-the-art tree-based methods have about log *N* times larger asymptotic access overhead. Yet, tree-based methods are still of practical interest because they tend to have much smaller access overhead constant factors than the hierarchical methods. The access overheads of tree-based methods also constrain the worst case while those of the hierarchical methods are often amortized. Although there are techniques to achieve competitive worst-case access overhead via the hierarchical methods [13, 19], they tend to be complex and add further constant factors to the performance bounds.

In the past, the ORAM research community has focused mainly on reducing the access overhead because it was the biggest obstacle to applying ORAM in practice. However, some recent studies have achieved practical access performances [16, 22] by combining tree-based ORAM with special hardware. For example, the PHAN-TOM secure processor system [16] supports access pattern-hiding SQL queries with a time overhead of 1.2–6 × compared to the standard insecure version. Thus, at least for tree-based ORAM, exploration of aspects other than access overhead is beginning to make sense.

3 The square root method [7] was the first non-trivial ORAM and is the origin of some of the ideas underlying the hierarchical methods.

We describe the recent development of techniques for reducing the number of physical cells *N* of the tree-based ORAM to (1 + *o*(1))*N* [17]. Note that there also is a space overhead originating from the cell size: typically, *B* = *B* + *c* lg *N* where *c* is a small constant such as 2. We ignore the cell size overhead and focus on the cell number because the typical value of *B* is 128 bytes in the secure processor setting and the overhead with respect to the cell size is just a few percent.

#### *9.2.3 Tree-Based Methods*

The tree-based method of Stefanov et al. [27] works as follows. The server organizes the physical cells into a complete binary tree with *N* leaves where each node is a bucket—a container that can accommodate a constant number of virtual cells. Each virtual cell has a position label—an integer in [*N*]—and a virtual cell with position label *i* is stored either in some bucket on the path from the root to the *i*th leaf or in a stash, which is a container in the user's local memory that can accommodate a small number of virtual cells. Let *vi* be the *i*th virtual cell and let *pi* be the position label of *vi* . Suppose the user maintains *pi* for all *i* ∈ [*N*] in local memory. This requires (*N*) user space but simplifies the exposition. We will reduce the user space later. To access *vi* , the user retrieves all of the blocks on the path from the root to the *pi* th leaf. Let this path be *P*. At this point, *vi* must be in the stash. The user copies the value of *vi* to somewhere in its local memory for a read query or changes it to some other value for a write query. Then, the user updates *pi* to a fresh random value in [*N*]. After that, the user scans the buckets on *P* from the leaf to the root, and for each bucket, moves cells in the stash to the bucket greedily while respecting the position labels and the bucket capacity. See Fig. 9.1 for an example.

Sometimes, some virtual cells in the stash cannot be moved back to the tree. For example, if all cells are assigned the same position label, only (log *N*) physical cells can be used to store the virtual cells and thus, most virtual cells must end up in the stash. (Of course, if *N* is large, such an event happens only extremely rarely.) Stefanov et al. proved that if the bucket size is at least 5, the number of cells left in the stash after processing a query is exponentially small. Thus, if the stash size is ω(log *N*), the stash overflows during processing a polynomial number of queries with only negligible probability.

To reduce the user space for position labels, the user outsources the position labels using the same method recursively. That is, each position label is a lg *N* bit integer and the table for position labels of all virtual cells can be thought of as a RAM storing the *N*lg *N* -bit concatenation of the integers. Thus, the original problem of hiding the access pattern to RAM consisting of *N* cells, each of *B* bits, is reduced to hiding the access pattern to RAM consisting of *N*lg *N* /*B* cells, each of *B* bits. If, say, *B* ≥ 2 lg *N*, which is a completely reasonable assumption for all reasonable *N*, 4 the problem size (cell number) decreases exponentially and reaches *O*(1) after

4 Recall that the typical value of *B* is 128 bytes in secure processor applications.

*O*(lg *N*) levels of recursion. At that point, the user can store the *O*(1) size RAM locally terminating the recursion.

The tree at the top level of recursion has size (*N*) and the tree at higher recursion levels decreases exponentially. The access overhead is proportional to the sum of the heights of the trees at all recursion levels, which is *O*(log2 *N*). The server space is proportional to the sum of the sizes of the trees at all recursion levels, which is (*N*).

Although each recursion level requires a stash, the numbers of cells left in those stashes are independent and it turns out that the total number of cells left in all stashes is still exponentially small. Thus, *f* (*N*) user space is enough for any *f* (*N*) ∈ ω(log *N*).

#### *9.2.4 Succinct Construction*

The constant factor hidden in the (*N*) server space bound of the method described above is about 10: the top-level tree has 2*N* nodes each of capacity 5 while the size of the recursive trees is negligible because typically, *B* is much larger than lg *N*. (Theoretically, we assume *B* = ω(lg *N*).) Though one can reduce this constant factor to some extent by decreasing the tree height while tuning the bucket size, it is not possible to achieve a factor ≤ 2 while maintaining a meaningful stash overflow probability, at least using the currently known analysis techniques. This method also leads to prohibitively large access overhead as the server space becomes close to 2*N*. We now describe a method for achieving (1 + *o*(1))*N* server space with a modest sacrifice in access overhead [17].

#### **Fig. 9.2** Large leaf layout

The idea is to modify the layout of the tree at the top recursion level so that the leaf number is *N*/ lg1.4 *N* and the leaf size is lg1.4 *N* + lg1.3 *N* (see Fig. 9.2). It is obvious that the tree size is (1 + 1/ lg0.1 *N*)*N* while access overhead remains *O*(log2 *N*). We now explain why the stash overflow probability remains small.

Let *Ni* be the number of cells with position label *i* for *i* ∈ [*N*/ lg1.4 *N*]. Let the load of a bucket be the number of cells stored in the bucket. At each moment in the lifetime of the scheme, *Ni* follows the binomial distribution with parameters *N* and *N*/ lg1.4 *N* for each *i*. The probability that this becomes larger than lg1.4 *N* + lg1.3 *N* is negligible. Thus, no leaf becomes full while processing a polynomial number of queries. Under this assumption, the distribution of the internal bucket loads is dominated by the distribution of the loads of the corresponding *N*/ lg1.4 *N* − 1 internal buckets in the standard *N* leaf layout scheme described above. This is so because the internal buckets in the large leaf layout do not need to store the cells that overflow from the leaves. Thus, assuming no leaf becomes full, the stash overflow probability of the large leaf case is negligible. The same is true even without this assumption because there is only a negligible probability that the assumed case does not occur.

We can reduce the *N*/ lg0.1 *N* extra term on the tree size even further by "the power of two choices". That is, we give two random position labels to each virtual cell. One is primary, which determines the path on which the cell can reside, while the other is secondary, which is a dummy needed to hide the access pattern. Now, *Ni* is the number of virtual cells with primary position label *i*. We maintain *Ni* for all *i* in a sub-ORAM in the same way we store position labels in recursive ORAM. To access a virtual cell *v*, we retrieve all cells on the path from the root to the *p*1th leaf and the path from the root to the *p*2th leaf. We choose two random labels *p* 1, *p* 2 and let *p* 1 (resp. *p* 2) be the new primary (resp. secondary) label of *v* if *Np* 1 < *Np* 2 . Otherwise, we exchange the role of *p* 1 and *p* 2. We then scan the paths specified by the old labels and greedily move back the cells as in the previous method. Here, *Ni* is not binomial but concentrated much more tightly around the mean due to the effect of the two choices. Thus, the "head space" for each leaf can be much smaller than lg1.3 *N*, leading to a smaller tree size.

By tuning the parameters, the first technique (large leaf layout) alone can achieve about (1 + (log *N B* + √ 1 log *N* ))*N* server space while the second technique (two choices) decreases it to (1 + (log *N B* + log log *N* log2 *N* ))*N*. 5

#### *9.2.5 Open Problem*

It is unknown whether the optimal *O*(log *N*) access overhead and (1 + *o*(1))*N* server space can be achieved at the same time. There are two natural approaches for answering this question affirmatively:

#### **9.3 Wear Leveling**

#### *9.3.1 Problem*

Consider the case where you have RAM with the limitation that each cell state can be updated at most a certain number of times. Once the number of updates has reached the limit, the cell dies and you can no longer update it. The utility of the RAM quickly degrades as the cells start to die because the total amount of information that can be stored decreases, and it becomes cumbersome to manage which cells are still alive. Thus, the number of times you can support updates before cells start to die is of primary interest. This number depends heavily on the case. In the best case where the updates are uniform among the cells, you can perform *nL* updates where *n* is the number of cells and *L* is the number of times each cell can be updated. In contrast, in the worst case where all updates fall onto a particular cell, you can perform updates only *L* times. Wear leveling is the problem of prolonging the memory lifetime as much as possible while keeping the associated overhead, if any, as small as possible.

The system community has been studying wear leveling for decades. Historically, flash memory was the main motivation for studies conducted from the late 1980s to the mid 2000s [1, 6, 15]. Today, the main motivation for wear leveling comes

5 These bounds include the cell size overhead that we ignored in the main explanation for brevity.

from phase change memory (PCM), which is an emerging next-generation memory technology that has many features, including low latency, energy efficiency, and non-volatility [5]. Each PCM cell supports only 108–109 updates, which means that cells can start dying within minutes or even seconds if no effort is made to perform wear leveling. PCM differs from flash memory in certain important respects, such as latency, access granularity, and in-place write capability, and thus requires a different wear leveling formalization than flash memory.

Most existing studies on wear leveling are conducted mainly from a practical point of view. Often, they do not have a formal problem statement or rigorous theoretical analyses. While this might not be a serious problem if the only thing that matters is the performance, some relatively recent studies have repeatedly emphasized the security aspects of wear leveling [21, 23, 24, 28, 29]. In particular, it is important to take into account the case of malicious users who actively try to reduce the memory lifetime. (Consider, for example, a computing outsourcing service.)

Below, we describe a recent theoretical study that constructed a problem formalization to capture the wear leveling for PCM explained above, and the corresponding solutions [18].

The formal problem statement is as follows. There are two parties: the user and server. The server has three resources: physical RAM, wear-free memory, and private randomness. The physical RAM is RAM that consists of *N B*-bit cells while the wearfree memory is RAM that consists of a small number of *B*-bit cells. The user provides the server with adversarially chosen read/write queries to virtual RAM—a RAM consisting of *n b*-bit cells—and the server must respond to these queries "correctly." That is, each request is of the form (type,*i*, *v*) where type ∈ {read, write}, *i* ∈ [*n*], and *v* ∈ {0, 1}*b* and, for (read,*i*, ∗) where "∗" means that the third component is arbitrary, the server must return the last value written to the *i*th virtual cell (the *v* in the last query of the form (write,*i*, *v*)). The server not only needs to return the correct responses but also needs to support as many write queries as possible with high probability without updating any physical cell more than *L* times where *L* is a parameter. We assume *L* = *n*δ for some constant δ > 0. Equivalently, we define δ := log*L n* and assume it is a constant. This assumption is reasonable even though, in reality, *L* and *n* are independent, because *L* is 108–109 and log*L n* is at most 2 or 3 for reasonable *n*.

The performance metric for wear leveling includes the physical memory size, the wear-free memory size, the number of write queries supported, and the access overhead, which is the number of physical RAM accesses needed for each virtual RAM access. We say a wear-leveling scheme is "optimal" if it satisfies the following conditions (asymptotic notations are in terms of *n* → ∞):

6 We ignore the cell size overhead for brevity. Security Refresh [23] described below does not have any cell size overhead (*B* = *B*) while the method of Onodera and Shibuya [18] described after that has a cell size of *B* = *B* + 2lg *n* + 1.

**Fig. 9.3** Movement of cells in an epoch of security refresh. All the numbers are binary. *n* = 4(= *N*), *r*0 ="10", and *r*1 ="11". Solid arrows mean cell swaps while dotted arrows mean skipped cell swaps

#### *9.3.2 Security Refresh*

The wear leveling scheme of Sewong et al. [23] is optimal if *L* = *N*δ for δ > 1 while it is non-optimal (in fact, "far from" optimal) for δ < 1 [18].

In this method, *n* virtual cells are stored in the *N* = *n* physical cells in permuted order. The method works in epochs. At each epoch, two random lg *n*-bit integers *r*0 and *r*1 are maintained. (We assume *n* is a power of two for brevity.) At the start of an epoch, for each *i* ∈ [*n*], the *i*th virtual cell *vi* is stored in the *i* ⊕ *r*0th physical cell *Vi*⊕*r*0 where ⊕ means bit-wise XOR. During the epoch, each *vi* is moved from *Vi*⊕*r*0 to *Vi*⊕*r*1 . Note that the virtual cell stored in the destination *Vi*⊕*r*1 of *vi* is *vi*⊕*r*1⊕*r*0 and its destination is *Vi*⊕*r*0 ; that is, *vi* and *vi*⊕*r*0⊕*r*1 swap their positions. This is done as follows. For every *t* write queries processed where *t* is a parameter, we perform a remap subroutine. At the *i*th remap subroutine call in an epoch, we check if *i* < *i* ⊕ *r*0 ⊕ *r*1. If so, *vi* still is in *Vi*⊕*r*0 and thus, we swap the contents of *Vi*⊕*r*0 and *Vi*⊕*r*1 . Otherwise, *vi* is already in *Vi*⊕*r*1 and we skip swapping. The epoch ends after the *n*th remap subroutine finishes. At that point, each cell *vi* is stored in *Vi*⊕*r*1 . We update *r*0 to *r*1, and *r*1 to a fresh random lg *n*-bit integer. Now every *vi* is in *Vi*⊕*r*0 as required for the epoch start, and we restart another epoch at this point. See Fig. 9.3 for an example. To access *vi* , we access *Vi*⊕*r*0 if *vi* was already remapped in the epoch. (We have already seen how to check this.) Otherwise, we access *Vi*⊕*r*1 .

The non-trivial part of the analysis is the proof of a high-probability guarantee of memory lifetime. We outline the key points. Fix a physical cell and let *Xi* be the number of times it is updated during the *i*th epoch. We need to place a bound on the probability that the sum of *Xi*s deviates from its expected value. To do this, it suffices to bound the deviation of the sum of odd-indexed variables *X*1, *X*3,... from its expected value and do the same for the sum of even-indexed variables *X*2, *X*4,... separately. This is helpful because each *Xi* is a random variable that depends on *r*0,*r*1 in the *i*th epoch (and the queries) and thus, the odd-indexed variables *X*1, *X*3,... are independent of each other and so are the even-indexed variables. Regardless of the queries, *Xi* is bounded by the number of write queries processed in an epoch *tn*. Although this suggests the use of the Hoeffding inequality, it turns out that it does not work for the case δ < 2, essentially because the condition *Xi* ≤ *tn* alone does not capture the fact that some cell being updated many, say, ≈ *tn*, times in an epoch negatively affects the number of times other cells are updated in the epoch. To derive the bound for the case 1 <δ< 2, bound the second moment of *Xi* and apply the Bernstein inequality [4].

If the user tries to keep on updating *vi* continuously, one of *Vi*⊕*r*0 and *Vi*⊕*r*1 is updated *tn*/2 = (*n*) times during the first epoch, and this physical cell dies if δ < 1. Thus, this method is not optimal for δ < 1.

#### *9.3.3 Construction for Small Write Limit Cases*

We now briefly describe a method for achieving optimality for the case δ < 1, that is, the memory is large [18]. The idea is to prepare spare cells and remap the frequently updated cells to free spare cells adaptively. (We maintain the write counts of cells by appending a counter to each cell.) We store pointers to the new locations in the old locations to trace the remapped cells. To keep the number of pointers to follow small, we connect pointers in a manner that is similar to the DFS of a complete *d*-ary tree with *dh* leaves where *d*, *h* are parameters (see Fig. 9.4). As we continue to process write queries, the data structure gradually degrades: the free spare cells become scarce and the trees become saturated. To reset the degradation, we perform a Security Refresh-style mapping. That is, we treat the structure in Fig. 9.4 as residing in another RAM *u* and maintain a global mapping—a gradually changing one-to-one map between the cells of *u* and the physical cells *V*1, *V*2,... . Once we have globally remapped a cell of *u* corresponding to a tree root, we reset the "DFS" starting from that cell. For example, if we globally remap *ui* in state (5) of Fig. 9.4, we free *un*+1, *un*+4, resetting DFS for the tree from *ui* to the root. Garbage such as *un*+3 are also reclaimed sooner or later when they are globally remapped. To access *vi* , the tree path traversal in *u* starting from *i* is simulated translating between addresses in *u* and addresses in *V*.

Although analysis of the bound on memory lifetime is cumbersome, the same idea as the analysis of Security Refresh applies. Indeed, the core argument is easier because the Hoeffding bound suffices.

#### *9.3.4 Open Problem*

The access overhead of the method for the small write limit case described above is about 1/δ. It is easy to obtain amortized 1 + *o*(1) and worst-case (*n*) access overhead if we allow relatively large wear-free memory, for example, *O*(*n* ) for an appropriate constant 0 << 1. It seems possible and practically relevant to

**Fig. 9.4** Example evolution of *u*. *d* = *h* = 2. *u*1,..., *un* are default locations while the rest are spare cells. Each panel shows the state just after the thick-bordered cell was allocated because the cell previously storing its content was updated and the write count reached the threshold

achieve amortized 1 + *o*(1) and worst-case *O*(1) access overhead in this setting. A theoretically more interesting challenge is to give negative results that justify the use of such large wear-free memory.

#### **9.4 Conclusion**

We reviewed two recent studies on ORAM and wear leveling that achieve succinct space usage. Though these objects have totally different motivations and are studied in different communities, there are some similarities between them. As we mentioned in the introduction, several other concepts with similar flavors are known, including initializable RAM, memory checking, locally decodable code, and huge random objects. There are probably many more such enhanced RAM instances yet to be found, and trying to find them can be an avenue for making progress in studies of data structures.

#### **References**

**Open Access** This chapter is licensed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license and indicate if changes were made.

The images or other third party material in this chapter are included in the chapter's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the chapter's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

# **Part IV Sublinear Modelling**

# **Chapter 10 Review of Sublinear Modeling in Probabilistic Graphical Models by Statistical Mechanical Informatics and Statistical Machine Learning Theory**

#### **Kazuyuki Tanaka**

**Abstract** We review sublinear modeling in probabilistic graphical models by statistical mechanical informatics and statistical machine learning theory. Our statistical mechanical informatics schemes are based on advanced mean-field methods including loopy belief propagations. This chapter explores how phase transitions appear in loopy belief propagations for prior probabilistic graphical models. The frameworks are mainly explained for loopy belief propagations in the Ising model which is one of the elementary versions of probabilistic graphical models. We also expand the schemes to quantum statistical machine learning theory. Our framework can provide us with sublinear modeling based on the momentum space renormalization group methods.

# **10.1 Introduction**

Statistical machine learning frameworks using **probabilistic graphical models** are useful for many applications, including information communication technologies [1–3], compressed sensing [4, 5] and neural information processing systems [6–10] in data-driven sciences.

Most probabilistic graphical models belong to the exponential family [11] and can be regarded as classical spin systems in **statistical mechanical informatics**[12–17]. However, it is well known that many applicable formulations in data sciences as well as computational sciences can be reduced to combinatorial problems with some constraint conditions which can be regarded as an **Ising Model** in statistical mechanical informatics [18, 19]. Moreover, much interest has focused on applying **quantum annealing** as a novel high-speed optimization technology to massive optimization problems [20–24].

K. Tanaka (B)

Graduate School of Information Sciences, Tohoku University, Sendai, Japan e-mail<EMAIL_ADDRESS>
N. Katoh et al. (eds.), *Sublinear Computation Paradigm*, https://doi.org/10.1007/978-981-16-4095-7\_10

#### **10.2 Statistical Machine Learning**

In statistical machine learning, most of the mathematical frameworks for machine learning are based on maximum likelihood frameworks [25, 26] from statistical mathematical sciences. The important points are how to assume the prior distribution and the data generative probability distribution and how to express the joint probability between the parameters and the data vector. In this section, we explore maximum likelihood frameworks in terms of model selection and parameter selection from a given data vector.

## *10.2.1 Bayesian Statistics and Maximization of Marginal Likelihood*

Let us consider a graph specified by nodes and edges,(*V*, *E*), where *V* is the set of all nodes*i* and *E* is the set of all edges{*i*, *j*}. State variables*si* and *di* are associated with

$$\text{Each node } i. \text{ The vectors } \mathbf{s} = \begin{pmatrix} s\_1 \\ s\_2 \\ \vdots \\ s\_{|V|} \end{pmatrix} \text{ and } \mathbf{d} = \begin{pmatrix} d\_1 \\ d\_2 \\ \vdots \\ d\_{|V|} \end{pmatrix} \text{ correspond to the parameters.}$$

and the data vector, respectively. The state spaces of *si* and *di* are given by and (−∞, +∞), respectively. Now ρ(*d*|*s*,β) and *P*(*s*|α) which correspond to the data generative and prior models, respectively, are assumed to be as follows:

$$\rho(\mathbf{d}|\mathbf{s},\beta) = \prod\_{i \in V} \sqrt{\frac{\beta}{2\pi}} \exp\left(-\frac{1}{2}\beta \left(d\_i - s\_i\right)^2\right),\tag{10.1}$$