ForestDSH: A Universal Hash Design for Discrete Probability Distributions

TL;DR

This paper introduces ForestDSH, a novel hash method for high-dimensional classification based on known probabilistic models, outperforming existing nearest neighbor search techniques in speed, especially in mass spectrometry applications.

Contribution

The paper presents a new distribution sensitive hashing scheme using decision forests, with a theoretical complexity analysis and practical superiority over current methods.

Findings

Faster search times than existing methods in theory and simulations.

Achieves an order of magnitude speedup in mass spectrometry spectral library search.

Complexity scales with N^{λ*(P)} where λ*(P) is analytically derived.

Abstract

In this paper, we consider the problem of classification of $M$ high dimensional queries $y^1,\cdots,y^M\in B^S$ to $N$ high dimensional classes $x^1,\cdots,x^N\in A^S$ where $A$ and $B$ are discrete alphabets and the probabilistic model that relates data to the classes $P(x,y)$ is known. This problem has applications in various fields including the database search problem in mass spectrometry. The problem is analogous to the nearest neighbor search problem, where the goal is to find the data point in a database that is the most similar to a query point. The state of the art method for solving an approximate version of the nearest neighbor search problem in high dimensions is locality sensitive hashing (LSH). LSH is based on designing hash functions that map near points to the same buckets with a probability higher than random (far) points. To solve our high dimensional classification problem, we introduce distribution sensitive hashes that map jointly generated pairs $(x,y)\sim P$ to the same bucket with probability higher than random pairs $x\sim P^A$ and $y\sim P^B$, where $P^A$ and $P^B$ are the marginal probability distributions of $P$. We design distribution sensitive hashes using a forest of decision trees and we show that the complexity of search grows with $O(N^{\lambda^*(P)})$ where $\lambda^*(P)$ is expressed in an analytical form. We further show that the proposed hashes perform faster than state of the art approximate nearest neighbor search methods for a range of probability distributions, in both theory and simulations. Finally, we apply our method to the spectral library search problem in mass spectrometry, and show that it is an order of magnitude faster than the state of the art methods.

Equations

$$Arg\underset{x\in X}{\max }\mathbb{P}(y\mid x),$$

$$\mathbb{P}(y\mid x)=\prod _{s=1}^{S}p(y_s\mid x_s).$$

$$\text{mbox}Ifd(x,y)\le R,\text{mbox}then{h}^{\mathcal{A}}(x)=h^{\mathcal{B}}(y)\text{mbox}withprobabilityatleast{p}_1.$$

$$\text{mbox}Ifd(x,y)\ge cR,\text{mbox}then{h}^{\mathcal{A}}(x)=h^{\mathcal{B}}(y)\text{mbox}withprobabilityatmost{p}_2.$$

$$Prob(H^{\mathcal{A}}(x)\cap H^{\mathcal{B}}(y)\ne \varnothing \mid (x,y)\sim \mathbb{P})$$

$$Prob(H^{\mathcal{A}}(x)\cap H^{\mathcal{B}}(y)\ne \varnothing \mid (x,y)\sim \mathbb{Q})$$

$$\mathbb{P}:{\mathcal{A}}^S\times {\mathcal{B}}^S\to [0,1],$$

$$p^{\mathcal{A}}(x_s)$$

$$p^{\mathcal{B}}(y_s)$$

$$q(x_s,y_s)$$

$$\alpha$$

$$\beta$$

$${\gamma }^{\mathcal{A}}$$

$${\gamma }^{\mathcal{B}}$$

$$Se{q}^{\mathcal{A}}(root)$$

$$Se{q}^{\mathcal{B}}(root)$$

$$Se{q}^{\mathcal{A}}(f(v,a,b))$$

$$Se{q}^{\mathcal{B}}(f(v,a,b))$$

$$H_z^{\mathcal{A}}(x)$$

$$H_z^{\mathcal{B}}(y)$$

$$\Phi (f(v,a_i,b_j))$$

$${\Psi }^{\mathcal{A}}(f(v,a_i,b_j))$$

$${\Psi }^{\mathcal{B}}(f(v,a_i,b_j))$$

$$\Psi (v)$$

$$\Phi (v)$$

$$\Psi (v)$$

$${\Psi }^{\mathcal{A}}(v)$$

$${\Psi }^{\mathcal{B}}(v)$$

$$\alpha (G)$$

$$\beta (G)$$

$${\gamma }^{\mathcal{A}}(G)$$

$${\gamma }^{\mathcal{B}}(G)$$

$$c_{tree}|V(G)|$$

$$+\left(\frac{c_{hash}N}{\alpha (G)}+\frac{c_{hash}M}{\alpha (G)}+\frac{c_{insertion}N{\gamma }^{\mathcal{A}}(G)}{\alpha (G)}+\frac{c_{insertion}M{\gamma }^{\mathcal{B}}(G)}{\alpha (G)}+\frac{c_{pos}MN\beta (G)}{\alpha (G)}\right)\log \frac{1}{1-TP},$$

$$|V(G)|$$

$$\frac{\alpha (G)}{\beta (G)}$$

$$\frac{\alpha (G)}{{\gamma }^{\mathcal{A}}(G)}$$

$$\frac{\alpha (G)}{{\gamma }^{\mathcal{B}}(G)}$$

$$\alpha (G)$$

$$c_{hash}(\#bands)|X|+c_{insertion}\sum _{z=1}^{\#bands}\sum _{x\in X}|H_z^{\mathcal{A}}(x)|$$

$$+c_{hash}(\#bands)|Y|+c_{insertion}\sum _{z=1}^{\#bands}\sum _{y\in Y}|H_z^{\mathcal{B}}(y)|,$$

$$\displaystyle\mathbb{E}\big{(}\sum_{z=1}^{\#bands}\sum_{x\in X}|H_{z}^{\mathcal{A}}(x)|\big{)}$$

$$=$$

Full text

1 Introduction

Consider the problem of classifying a large number of high dimensional data $Y=\{y^{1},\cdots,y^{M}\}\subset\mathcal{B}^{S}$ into high dimensional classes $X=\{x^{1},\cdots,x^{N}\}\subset\mathcal{A}^{S}$, given a known joint probability distribution $\mathbb{P}(x,y)$, where $\mathcal{A}$ and $\mathcal{B}$ are discrete alphabets. Given a point $y\in Y$, the goal is to find the class $x\in X$ that maximize $\mathbb{P}(y\mid x)$;

[TABLE]

where $\mathbb{P}(y\mid x)$ is factorizable to i.i.d. components, i.e.,111Here, we assume that $\mathbb{P}(x)$ is also factorizable to i.i.d. components.

[TABLE]

In this paper, we refer to $X=\{x^{1},\cdots,x^{N}\}$ as data base points and $Y=\{y^{1},\cdots,y^{M}\}$ as queries. This problem has application in various fields, including the clustering of spectra generated by mass spectrometry instruments [1]. Consider the problem where there are billions of data points (mass spectra) and given a query spectrum $y$ the goal is to find the spectrum $x$ that maximize known probability distribution $\mathbb{P}(y\mid x)$ [2, 3]. This is similar to the nearest neighbor search problem [4, 5], where instead of minimizing a distance, we maximize a probability distribution, and has applications in machine learning [6], database querying [7] and computer vision [8, 9] .

A naive approach to solve this problem is to compute $\mathbb{P}(y\mid x)$ for each $x\in X$, and find the maximum. The run time for this algorithm is $O(NS)$ which is very slow when the number of classes, $N$, is large.

In order to address this problem where the number of classes is massive, multiclass classification methods have been established [10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22]. For example, in [16] an approach to construct a tree with logarithmic (in terms of number of classes) depth is presented and logarithmic train and test time complexity is achieved. However, these methods are limited to low-dimensional data (when the number of dimensions is less than $10$), and their complexity grows exponentially with the dimension. The methods in [10, 11, 12] speed up the prediction by rejecting a significant portion of classes for each query. However, all these methods suffer from the curse of dimensionality, i.e., when the dimension of the data increases, the complexity of these methods increases linearly with the number of classes.222The curse of dimensionality holds for non-deterministic probability distributions. When $p(y\mid x)$ is deterministic, i.e., when it takes only a single value with probability one, there is no curse of dimensionality. In this paper, we are interested in the case of non-deterministic probability distributions. Dimension reduction can also be utilized before applying nearest neighbor search algorithms to high dimensional data [23, 24, 25, 26]. Moreover, Structure Preserving Embedding (SPE) which is a low-dimensional embedding algorithm for Euclidean space is also used as a dimension reduction tool [27]. Set similarity in nearest neighbor search literature is studied extensively [28, 29]. In [28], the authors develop a data structure that solves the problem of approximate set similarity search under Braun-Blanquet similarity $B(x,y)=\frac{|x\cap y|}{\max(|x|,|y|)}$ in sub-quadratic query time.

A similar problem has been investigated in the field of nearest neighbor search. In this problem, given a set of points in a database, the goal is to find the point in the database that is closest to a query point. A popular approach to this problem is locality sensitive hashing [30, 31]. This approach solves $\epsilon$-approximate nearest neighbor search problem by designing a family of hashes in a way that near points are hashed to the same bucket with probability much higher than random points. In $\epsilon$-approximate nearest neighbor search problem, given a query point $y$, the goal is to find $x\in X$ for which $d(x,y)\leq(1+\epsilon)d(x^{\prime},y)$ for all $x^{\prime}\in X$ and $X$ is the set of all feasible points [30, 32]. One of the most popular methods for approximate nearest neighbor search is LSH [30, 31, 33, 34, 35, 28]. For any metric space $\mathcal{M}=(M,d)$, a family of functions $h:{\mathcal{M}}\to S$ is $\mathcal{F}(R,cR,{p}_{1},{p}_{2})$-sensitive if for any two points $x,y\in{\mathcal{M}}$:

[TABLE]

A family is interesting when $p_{1}>p_{2}$. In the case of hamming distance, i.e., when data are in the form of $d$-dimensional vectors from ${\{0,1\}}^{d}$, the family of hashes may be defined simply as $H=\cup_{i=1}^{d}\{h:\{0,1\}^{d}\to\{0,1\}\mid h(x)=x_{i}\}$ where $x_{i}$ is $i$-th coordinate of $x$. From (3) and (4), we conclude that $p_{1}=1-\frac{R}{d}$ and $p_{2}=1-\frac{cR}{d}$333Note that $d(x,y)\leq R$ is equivalent to $x$ and $y$ are different in at most $R$ coordinates..

In this paper, in addition to going from minimizing a distance to maximizing a probabilistic measure, our approach differs from the classic LSH in the following way. The proposed hashes in this paper are defined to be a subset of integers instead of an integer and our goal is to ensure that444In fact, to control the complexity, two more conditions are defined in Definition 2. ${\mathbb{P}(x,y)}$ is joint probability distribution while ${\mathbb{P}^{\mathcal{A}}(x)}$ and ${\mathbb{P}^{\mathcal{B}}(y)}$ are marginal probability distributions of ${\mathbb{P}(x,y)}$. ${\mathbb{Q}(x,y)}$ is also defined as ${\mathbb{Q}(x,y)}={\mathbb{P}^{\mathcal{A}}(x)}{\mathbb{P}^{\mathcal{B}}(y)}$.

[TABLE]

In words, these hashes hash the jointly-generated pairs of points to the same buckets with probabilities higher than $\alpha$ while hash random pairs of points to the same buckets with probabilities lower than $\beta$. Note that, in this case for data points $x$ and $y$, collision happens when we have $H^{\mathcal{A}}(x)\cap H^{\mathcal{B}}(y)\neq\varnothing$, while in the classic LSH, $x$ and $y$ have collision if $H^{\mathcal{A}}(x)=H^{\mathcal{B}}(y)$. The idea of defining hash collision as the intersection of hash sets has been previously proposed in [28].

Currently, the locality sensitive hashing approach does not generalize to the cases where the triangle inequality, i.e., $d(x,y)\leq d(x,z)+d(z,y)$ does not hold [29, 30, 31, 36, 37, 38, 39, 40]. These papers are based on balls on some specific points in the metric space and the notion of balls is well defined only for metric spaces that satisfy triangle inequality. Recently, high dimensional approximate nearest neighbors in a known probabilistic distribution setting have been investigated [32, 41]. However, currently it is not possible to design efficient algorithms based on these methods, due to the large number of parameters involved.

The problem of finding high dimensional approximate nearest neighbors in a known probabilistic setting using a bucketing tree algorithm has been studied previously in [41, 42]. Dubiner uses a strategy to hash the data points from an arbitrary joint probability distribution into the leafs of the tree in a way that the paired data collide with a probability higher than the random pairs. Here, paired data refers to pairs coming from joint probability distribution, i.e., $(x,y)\sim\mathbb{P}$ random pairs are the pairs coming from an independent probability distribution, i.e., $x\sim{\mathbb{P}}^{\mathcal{A}}$ and $y\sim{\mathbb{P}}^{\mathcal{B}}$. However, the algorithm introduced in Dubiner requires solving computationally intractable optimizations, making it impossible to implement (e.g. see equation (126) from [41]). However, in this paper, for the specific case where the distribution for any $s\in\{1,2,\cdots,S\}$ is $p(x_{s}=0,y_{s}=0)=p(x_{s}=1,y_{s}=1)=\frac{p}{2},p(x_{s}=1,y_{s}=0)=p(x_{s}=1,y_{s}=0)=\frac{1-p}{2}$, our theoretical and practical results are compared to [41] for the range of $0\leq p\leq 1$.

In this paper, we propose to solve (1) by defining a family of distribution sensitive hashes satisfying the following property. They hash the jointly-generated pairs of points to the same buckets with probabilities much higher than random pairs. We further design an algorithm to solve (1) in sub-linear time using these families of hashes. Next, a method to find optimal family of hashes is presented to achieve minimum search complexity using multiple decision trees. Note that, these decision trees have the same tree structure while they apply to different permutations $\{1,2,\cdots,S\}\rightarrow\{1,2,\cdots,S\}$ of the data. This way, we design forest of decision trees where each decision tree captures a very small ratio $\alpha$ of true pairs for some $\alpha\in\mathbb{R}^{+}$ and by recruiting $\#bands=O(\frac{1}{\alpha})$ independently permuted decision trees we can reach near perfect recovery of all the true pairs. In this paper, we refer to each decision tree as a band and $\#bands$ is referred to as the number of bands.

The main idea is that we construct decision-tree hashes, in a way that the chance of true pairs hashing to the same leaf nodes in these decision trees is higher than random points (Figure $1$). The decision tree is built in a way that the ratio $\frac{\mathbb{P}(x,y)}{\mathbb{Q}(x,y)}$ is higher than a minimum threshold, while the ratios $\frac{\mathbb{P}(x,y)}{\mathbb{P}^{\mathcal{A}}(x)}$ and $\frac{\mathbb{P}(x,y)}{\mathbb{P}^{\mathcal{B}}(y)}$ and the number of nodes in the graph are lower than a maximum thresholds, see Algorithm 4. We further determined the optimal tree among many trees that can be constructed in this way. Two theorems are presented here on the complexity of the decision tree built in Algorithm 4.

No decision tree exists with the overall search complexity below $O(N^{\lambda^{*}})$ where $\lambda^{*}$ is derived analytically from the probability distribution $\mathbb{P}(x,y)$. 2. 2.

The decision tree construction of Algorithm 4 described in (57) results in the overall search with complexity $O(N^{\lambda^{*}})$.

Our results show that our approach, Forest-wise distribution sensitive hashing (ForestDSH), provides a universal hash design for arbitrary discrete joint probability distributions, outperforming the existing state of the art approaches in specific range of distributions, in theory and practice. Moreover, we applied this method to the problem of clustering spectra generated from mass spectrometry instruments.

An alternative strategy for solving (2) is to reformulate the problem as minimum inner product search problem (MIPS) [35] by transferring the data points into a new space. However, as we show in Section 6 the transferred data points are nearly orthogonal to each other making it very slow to find maximum inner product using the existing method [35].

Note that, the algorithms presented in this paper are based on the assumption that the true pairs are generated from a known distribution $\mathbb{P}$. The distribution $\mathbb{P}$ can be learned from a training dataset of true pairs. In practice, training data can be collected by running a brute force search on smaller data sets or portion of the whole data. For example, in case of mass spectrometry search, we collect training data by running brute-force search on small portion of our data [1]. Note that, we can never perfectly learn the probability distribution that the data is generated from. In Theorem 3, we prove that our algorithm is robust to noise and our experiment confirms this.

Notation: The cardinality of a set $A$ is denoted as $|A|$. The sets $\mathbb{N}$ and $\mathbb{R}$ stand for the sets of natural and real numbers, respectively. We use $\mathbb{P}(\cdot)$ and $\mathbb{Q}(\cdot)$ to denote the probability function $\mbox{Prob}(\cdot)$. $\mathbb{E}(\nu)$ denotes the expected value of the random variable $v$. Moreover, we use the notation $f(x)=O(g(x))$, if $\limsup_{x\rightarrow\infty}\frac{|f(x)|}{g(x)}<\infty$ and the notation $f(x)=\Omega(g(x))$, if $\limsup_{x\rightarrow\infty}\frac{|f(x)|}{g(x)}>0$. In this paper, $\log x$ is computed to the base $e$.

2 Definitions

Definition 1

Define the $S$-dimensional joint probability distribution $\mathbb{P}$ as follows.

[TABLE]

where $\mathcal{A}=\{a_{1},a_{2},\cdots,a_{k}\}$, $\mathcal{B}=\{b_{1},b_{2},\cdots,b_{l}\}$, $x=(x_{1},x_{2},\cdots,x_{S})\in\mathcal{A}^{S}$ and $y=(y_{1},y_{2},\cdots,y_{S})\in\mathcal{B}^{S}$ where $k,l<\infty$. On the other hand, we assume that the probability distribution function ${p}(a,b)$ is independent of $s$ and satisfies $\sum_{a\in[k],b\in[l]}p_{a,b}=1$. Similarly, we define the marginal probability distributions $\mathbb{P}^{\mathcal{A}}:\mathcal{A}^{S}\rightarrow[0,1]$, $\mathbb{P}^{\mathcal{B}}:\mathcal{B}^{S}\rightarrow[0,1]$ and $\mathbb{Q}:\mathcal{A}^{S}\times\mathcal{B}^{S}\rightarrow[0,1]$ as $\mathbb{P}^{\mathcal{A}}(x)=\prod_{s=1}^{S}{p}^{\mathcal{A}}(x_{s})$, $\mathbb{P}^{\mathcal{B}}(y)=\prod_{s=1}^{S}{p}^{\mathcal{B}}(y_{s})$ and $\mathbb{Q}(x,y)=\prod_{s=1}^{S}{q}(x_{s},y_{s})$, where

[TABLE]

We use $p_{ij}$ instead of $p(a_{i},b_{j})$ for simplicity. Moreover, $p_{i}^{\mathcal{A}}$, $p_{j}^{\mathcal{B}}$ and $q_{ij}$ are defined as $\sum_{j=1}^{l}p_{ij}$, $\sum_{i=1}^{k}p_{ij}$ and $q(a_{i},b_{j})$, respectively. Finally, we use compact notations $P=[p_{ij}]$ and $Q=[q_{ij}]$ as the $k\times l$ matrices with $p_{ij}$ and $q_{ij}$ in their $i$-th row and $j$-th column, respectively.

We define family of distribution sensitive hashes as follows.

Definition 2

[Family of Distribution Sensitive Hashes] Assume that the four parameters $\alpha$, $\beta$, $\gamma^{\mathcal{A}}$ and $\gamma^{\mathcal{B}}$ along with the probability distributions $\mathbb{P}$, $\mathbb{P}^{\mathcal{A}}$, $\mathbb{P}^{\mathcal{B}}$, $\mathbb{Q}$ and finite set $V_{Buckets}$ are given where $\mathbb{P}^{\mathcal{A}}$ and $\mathbb{P}^{\mathcal{B}}$ are marginal distributions of $\mathbb{P}$ and $\mathbb{Q}=\mathbb{P}^{\mathcal{A}}\mathbb{P}^{\mathcal{B}}$. A family of hashes $H_{z}^{\mathcal{A}}:\mathcal{A}^{S}\rightarrow 2^{V_{Buckets}}$ and $H_{z}^{\mathcal{B}}:\mathcal{B}^{S}\rightarrow 2^{V_{Buckets}}$ is called $(\mathbb{P},\alpha,\beta,\gamma^{\mathcal{A}},\gamma^{\mathcal{B}})$-distribution sensitive, if the following hold

[TABLE]

where $1\leq z\leq\#bands$ and $\#bands$ represent the number of bands, while $V_{Buckets}$ represent a set of indices for the buckets. We show how to choose $\#bands$ in (51) in Appendix C and how to select $V_{Buckets}$ in Algorithm 3. Intuitively, $\alpha$ represents the chance of a true pair falling in the same bucket, while $\beta$ represents the chance of random pairs falling in the same bucket. We will show that $\gamma^{\mathcal{A}}$ and $\gamma^{\mathcal{B}}$ represent the complexity of computing which buckets the data points fall into. In the next section, we will describe how the families of distribution sensitive hashes can be used to design an efficient solution for (1). Note that, in [28] definitions similar to Definition 2 have been made in order to speed up set similarity search. How the data points are mapped through a decision tree is sketched in Figure 1. The decision tree and the set of buckets are explained in Algorithm 4. Finally, in Figure 2, we show how the rates of positive call for true and random pairs, i.e., $\alpha$ and $\beta$ are derived. The average number of times that data points and queries fall in buckets, i.e., $\gamma^{\mathcal{A}}$ and $\gamma^{\mathcal{B}}$ are also derived and shown in this figure.

Remark 1

In classic LSH, we have $|H^{\mathcal{A}}(x)|=1$ and $|H^{\mathcal{B}}(y)|=1$. Therefore, $\gamma^{\mathcal{A}}=\gamma^{\mathcal{B}}=1$. Our approach is more flexible in the following two ways. First, we allow for $\gamma^{\mathcal{A}}$ and $\gamma^{\mathcal{B}}$ to be larger than one. Moreover, even in the case of $\gamma^{\mathcal{A}}=\gamma^{\mathcal{B}}=1$, our method can optimize over the larger set of possible hashes.

Here we present a general algorithm on how ForestDSH algorithm works.

3 Forest Distribution Sensitive Hashing Algorithm

In this section, we assume an oracle has given us a family of distribution sensitive hashes $H_{z}^{\mathcal{A}},H_{z}^{\mathcal{B}},\forall z:1\leq z\leq\#bands$ that satisfies $(\eqref{family1}-\eqref{family5})$. Inspired by LSH method, we present Algorithm 2 for solving (1) using this family.

Remark 2

Note that, the number of times $\mathbb{P}(y\mid x)$ is computed in the brute force method to solve (1) is $|X|$. Note that, in the optimization problem (1), the point $y\in Y$ is given. The goal of Algorithm 2 is to solve (1) with a much smaller number of comparisons than the brute force.

Remark 3

In the special case when the hashes are coming from a decision tree, we analyze the complexity of Algorithm 2 in Section 5. We show that the number of positive calls in Algorithm 2 is proportional to $\beta$, while the complexity of computing $|H_{z}^{\mathcal{A}}(x)|$ and $|H_{z}^{\mathcal{B}}(y)|$ are proportional to $\gamma^{\mathcal{A}}$ and $\gamma^{\mathcal{B}}$. Moreover, the chance of true pairs being called positive grows with $\alpha$. Therefore, in the next sections, we attempt to design buckets such that $\alpha$ is maximized, while $\beta$, $\gamma^{\mathcal{A}}$ and $\gamma^{\mathcal{B}}$ are minimized.

Now, the question is how we can design these families in a way to minimize the complexity, and efficiently map data points to these families. We investigate these questions in Sections 4 and 4.1.

4 Designing ForestDSH Using a Decision Tree Structure

In the previous section, we assumed that an oracle has given us a family of distribution sensitive hashes. In this section, we design buckets that satisfy $(\eqref{family1}-\eqref{family5})$ using a forest of decision trees with the same structure. Here, we focus on the probability distributions that can be factorized as the product of i.i.d. components.

Each of our decision trees recovers ratio $\alpha$ of true pairs and by recruiting $\#bands=O(\frac{1}{\alpha})$ decision trees we can recover nearly all true pairs. This can be more efficient than using a single decision tree classifier as achieving near perfect true pair recovery by the single decision tree would require near brute-force complexity. By allowing $\alpha<1$, we can select decision tree that avoid paths with low $\mathbb{P}(x,y)$ resulting in complexities much lower than the brute-force search. Note that $\alpha$ is the true positive rate, e.g. the probability that a true pair $(x,y)$ fall in the same bucket, therefore we always have $\alpha\leq 1$, see Figure 2.

Assume that a decision tree $G=(V,E,f)$ is given where $V$ is the set of nodes, $E$ is the set of edges, $V_{l}\subset V$ is the set of leaf nodes in the decision tree and $f:V/V_{l}\times\mathcal{A}\times\mathcal{B}\rightarrow V$ is the decision function, see Figure 1. For the two nodes $v_{1},v_{2}\in V$, $v_{1}$ is called an ancestor of $v_{2}$ if $v_{1}$ falls within the path from $v_{2}$ to root. In this case, $v_{2}$ is called a descendant of $v_{1}$. Furthermore, assume that a subset of leaf nodes $V_{Buckets}\subset V_{l}$ is given and at depth $s$ in the decision tree, the decisions depend only on $x_{s}$ and $y_{s}$ where $x=(x_{1},\cdots,x_{S})$ and $y=(y_{1},\cdots,y_{S})$ are $S$-dimensional data point and query, respectively. We define functions $Seq^{\mathcal{A}}:V\rightarrow\cup_{s=0}^{S}\mathcal{A}^{s}$ and $Seq^{\mathcal{B}}:V\rightarrow\cup_{s=0}^{S}\mathcal{B}^{s}$ recursively as:

[TABLE]

where $[S,a]$ stands for the concatenation of string $S$ with character $a$, i.e., for a string $S=s_{1},\cdots,s_{n}$ of length $n$, $[S,a]$ would be $s_{1},\cdots,s_{n},a$ which is a string of length $n+1$. Moreover, given permutations $p_{z}:\{1,\cdots,S\}\rightarrow\{1,\cdots,S\},1\leq z\leq\#bands$, $x=(x_{1},\cdots,x_{S})$ and $y=(y_{1},\cdots,y_{S})$ define $perm_{z}(x)=(x_{p_{z}(1)},x_{p_{z}(2)},\cdots,x_{p_{z}(S)})$ and $perm_{z}(y)=(y_{p_{z}(1)},y_{p_{z}(2)},\cdots,y_{p_{z}(S)})$. Finally, the family of buckets $H_{z}^{\mathcal{A}}(x)$ and $H_{z}^{\mathcal{B}}(y)$ are defined as

[TABLE]

Note that the permutation $perm_{z}$ is a deterministic function of $z$, i.e., at each band the same permutation is used to randomly permute the data points $x$ and $y$. These permutations are chosen before mapping the data points. In Figure 3, we show how both $x$ and $y$ data points are first permuted randomly (using the same permutation) and then are mapped to the buckets in decision trees. Note that, here we call a pair positive, if they fall into the same bucket in at least one of the bands. Now, we show that these hashes are distribution sensitive.

Definition 3

The functions $\Phi:V\rightarrow\mathbb{R}$, ${\Psi}^{\mathcal{A}}:V\rightarrow\mathbb{R}$ and ${\Psi}^{\mathcal{B}}:V\rightarrow\mathbb{R}$ are defined as follows. At root, $\Phi({root})=1$, ${\Psi}^{\mathcal{A}}({root})=1$ and ${\Psi}^{\mathcal{B}}({root})=1$, and for $a_{i}\in\mathcal{A}$, $b_{j}\in\mathcal{B}$ and $v\in V$, $\Phi(v)$, ${\Psi}^{\mathcal{A}}(v)$, and ${\Psi}^{\mathcal{B}}(v)$ are defined recursively as

[TABLE]

Moreover, $\Psi:V\rightarrow\mathbb{R}$ is defined as

[TABLE]

Lemma 1

The following properties hold:

[TABLE]

Remark 4

Note that the left side of $(\eqref{fam1}-\eqref{fam2})$ is independent of $z$ while it seems like the right side depend on band $z$. The proof of Lemma 1 in Appendix A shows that in fact it is independent of band $z$.

Lemma 2

For any decision tree $G$, satisfying the condition that for any pair of buckets $v_{1},v_{2}\in V_{Buckets}$, $v_{1}$ is not an ancestor or descendant of $v_{2}$, $H_{z}^{\mathcal{A}}(x)$ and $H_{z}^{\mathcal{B}}(y)$ defined in (19) and (20) are $(\mathbb{P},\alpha(G),\beta(G),\gamma^{\mathcal{A}}(G),\gamma^{\mathcal{B}}(G))$-sensitive where

[TABLE]

Proofs of Lemmas 1 and 2 are relegated to Appendix A. So far, we showed how to design ForestDSH using a decision tree structure. However, it is not yet clear how to map data points to these buckets. In the next section, we investigate this and provide an algorithm to design optimal decision trees.

4.1 Mapping Data Points

In the previous sections, we presented Algorithm 2 for solving (1) where we need to compute $H_{z}^{\mathcal{A}}(x)$ and $H_{z}^{\mathcal{B}}(y)$ by mapping data points to the buckets. We did not clarify how this mapping can be done efficiently. We present Algorithm 3 for mapping data points to the buckets using hash-table search, see Figure 2.

Remark 5

Queries are mapped to buckets using hash-table search similar to data points through Algorithm 3.

Note that we slightly modify the hash-table to search for values that are prefix of a query, rather than being exactly identical. In Section 5, we show that the complexity of Algorithms 2 and 3 can be formulated as:

[TABLE]

where $c_{tree}$, $c_{hash}$, $c_{insertion}$ and $c_{pos}$ are constants not depending on $N$. For the intuition behind (33), note that the first term $c_{tree}|V(G)|$ stands for the time required for calculating and storing the tree. The second and third terms, i.e., $\left(\frac{c_{hash}N}{\alpha(G)}+\frac{c_{hash}M}{\alpha(G)}\right){\log\frac{1}{1-TP}}$ denote the time needed for inserting data points to the hash-table. The fourth and fifth terms, i.e., $\left(\frac{c_{insertion}N\gamma^{\mathcal{A}}(G)}{\alpha(G)}+\frac{c_{insertion}M\gamma^{\mathcal{B}}(G)}{\alpha(G)}\right){\log\frac{1}{1-TP}}$ stand for the time required for mapping the data points from the hash-table to buckets. Finally, the last term, i.e., $\left(\frac{c_{pos}MN\beta(G)}{\alpha(G)}\right){\log\frac{1}{1-TP}}$ is the time of brute-force checking within each bucket.

In order to bound (33) with $O(N^{\lambda})$ for some $\lambda\in\mathbb{R}^{+}$, it is necessary and sufficient to find a tree $G$ that satisfies the following constraints:

[TABLE]

where ${\delta}=\frac{\log M}{\log N}$.

Theorem 1

Complexity of Algorithms 2 and 3 is equal to (33). Moreover, there exists a decision tree $G$ that is optimal and satisfies $(\eqref{order1}-\eqref{order3})$ for $\lambda$ defined in Definition $4$.

Theorem 1 is proved in three steps: $(a)$. In Section 5, we prove that the complexity is equal to (33). $(b)$. Theorem 2 shows that the tree $G$ constructed by Algorithm 4 and presented in Section 6 satisfies ((34)-(38)) for $\lambda$ defined in Definition 4. $(c)$. Theorem 3 shows that this is optimal.

5 Complexity Analysis

In this section, we bound the complexity of Algorithms 2 and 3 by (33). Note that, the complexity of Algorithms 2 and 3 is the summation of the following terms.

Tree construction complexity. $c_{tree}|V(G)|$ is the tree construction complexity where $c_{tree}$ is a constant representing per node complexity of constructing a node and $|V(G)|$ is the number of nodes in the tree. 2. 2.

Data mapping complexity. The complexity of this hash-table search grows with

[TABLE]

where $c_{hash}$ and $c_{insertion}$ represent complexity of insertion in the hash-table and insertion in the buckets, respectively. 3. 3.

Complexity of checking positive calls. (\#bands)c_{pos}\sum_{x\in X,y\in Y}|H_{z}^{\mathcal{A}}(x)\cap H_{z}^{\mathcal{B}}(y)|\big{)} is the complexity of checking positive calls where $c_{pos}$ is a constant representing the complexity of computing $\mathbb{P}(y\mid x)$ for a positive. Note that, $c_{tree}$, $c_{hash}$, $c_{insertion}$ and $c_{pos}$ are constants not depending on $N$. From $(\eqref{family1}-\eqref{family5})$, we have

[TABLE]

Note that, the total number of collision for random pairs is the sum of number of buckets that they intersect at. Therefore, we conclude that

[TABLE]

where (45) is concluded as ${\mid H^{\mathcal{A}}(x)\cap H^{\mathcal{B}}(y)\mid}\leq 1$.

Now, the question is how we can select $\#bands$ such that the true positive rate, defined as the ratio of true pairs that are called positive is high. In each band, the chance of a pair $(x,y)\sim\mathbb{P}$ being called positive is computed as $\alpha(G)=\sum_{v\in V_{Buckets}(G)}\Phi(v)$. Therefore, the overall true positive rate can be computed as:

[TABLE]

Using (50), and the inequality $(1-x)^{\frac{c}{x}}<e^{-c}$, the minimum possible value of $\#bands$ to ensure true positive rate $TP$ can be computed as

[TABLE]

where $\lceil r\rceil$ stands for the smallest integer greater than or equal to $r$. Therefore, the total complexity is equal to (33).

6 Constructing optimal decision trees for ForestDSH

In this section, we present an algorithm to design decision trees with complexity $O(N^{\lambda^{*}})$, where $\lambda^{*}$ is defined below, and we show that it is the optimal decision tree.

Definition 4

Given probability distributions $P=[p_{ij}]$ and $Q=[q_{ij}]$, $1\leq i\leq k,1\leq j\leq l$, and number of queries and classes $M$ and $N$ define $\delta=\frac{\log M}{\log N}$ and

[TABLE]

Remark 6

*For any probability distribution $\mathbb{P}$, the parameters $\mu^{*}$, $\nu^{*}$, $\eta^{*}$ and $\lambda^{*}$ can be derived numerically from Algorithm 5 in  Appendix  B. The intuition  behind  the  definition  of   $\mathcal{I}$ and $(\mu^{*},\nu^{*},\eta^{*})$ is that in Lemma 4 in Section 6.1 we show that for any decision tree $G$ and the variables $\Phi(v)$, ${\Psi}^{\mathcal{A}}(v)$, ${\Psi}^{\mathcal{B}}(v)$ and $\Psi(v)$ defined in (21)-(24), we have

\sum_{v\in V_{Buckets}(G)}{\big{(}\Phi(v)\big{)}}^{1+\mu+\nu-\eta}{\big{(}{\Psi}^{\mathcal{A}}(v)\big{)}}^{-\mu+\eta}{\big{(}{\Psi}^{\mathcal{B}}(v)\big{)}}^{-\nu+\eta}{\big{(}\Psi(v)\big{)}}^{-\eta}\leq 1 if $(\mu,\nu,\eta)\in\mathcal{I}$. Moreover, in proof of Theorem 2 we show that $(\mu^{*},\nu^{*},\eta^{*})$ are Lagrangian multipliers in an optimization problem to minimize the search complexity in Algorithm 2 while retaining a nearly perfect recovery. Consider the optimal decision tree and all the nodes $v$ satisfying*

Depth of the node $v$ is $n^{*}$. 2. 2.

For any $1\leq i\leq k,1\leq j\leq l$, the ratio of times we have $a_{i}$ at $s$-th position of $x$ and $b_{j}$ at $s$-th position of $y$ in all $1\leq s\leq n^{*}$ is $r_{ij}$.

The node $v$ or one of its ancestors is designated as a bucket by Algorithm 4. This is proved in Section 6.2.

In Algorithm 4, we provide an approach for designing decision trees with complexity $O(N^{\lambda^{*}})$. The algorithm starts with the root, and at each step, it either accepts a node as a bucket, prunes a node, or branches a node into $kl$ children based on the following constraints555If for any node the constraints for accepting as a bucket and pruning hold simultaneously, the algorithm accepts the node as a bucket. :

[TABLE]

Note that, for the cases where $p_{ij}=0$, we simply remove that branch, therefore, we assume all the $p_{ij}$ are positive real numbers. In Section 5, we prove that in order to bound the complexity in (33) with $O(N^{\lambda})$ for some $\lambda\in\mathbb{R}^{+}$, it is necessary and sufficient to find a decision tree $G$ that satisfies the constraints $(\eqref{order1}-\eqref{order3})$.

Remark 7

In Theorem 3, we prove that the decision tree construction of Algorithm 4 results in a tree with complexity $O(N^{\lambda^{*}})$ by setting $C_{1}=C_{2}=C_{3}=p_{0}q_{0}$ where $p_{0}$ and $q_{0}$ are defined as $\prod_{i,j}p_{ij}$ and $\min(\prod_{i,j}q_{ij},\prod_{i}{(p_{i}^{\mathcal{A}})}^{l},\prod_{j}{(p_{j}^{\mathcal{B}})}^{k})$. Pessimistic lower bounds $C_{1}=C_{2}=C_{3}=p_{0}q_{0}$ are derived in proof of Theorem 3 and in practice $C_{1}$, $C_{2}$ and $C_{3}$ are chosen in a way that the best complexity is achieved.

Example 1

Here, we focus on the case where $\mathcal{A}=\{0,1\},\mathcal{B}=\{0,1\}$, $P=\begin{bmatrix}0.4&0.3\\ 0.1&0.2\end{bmatrix}$, $Q=\begin{bmatrix}0.35&0.35\\ 0.15&0.15\end{bmatrix}$, $\delta=1$, and $M=N=4$. From Algorithm 5, we have $\mu^{*}=12.0791$, $\nu^{*}=13.4206$, $\eta^{*}=11.0959$, $\lambda^{*}=1.7203$. The decision tree is constructed from Algorithm 4 and is depicted in Figure 4. The nodes in the tree that are selected as bucket, i.e., satisfying $\frac{\Phi(v)}{\Psi(v)}\geq C_{1}{N^{1+\delta-\lambda^{*}}}$, are shown with a green check mark, and the nodes pruned out, satisfying either $\frac{\Phi(v)}{{\Psi}^{\mathcal{A}}(v)}\leq C_{2}N^{1-\lambda^{*}}$ or $\frac{\Phi(v)}{{\Psi}^{\mathcal{B}}(v)}\leq C_{3}N^{\delta-\lambda^{*}}$ are shown with a red cross. For the non-leaf (intermediate) nodes, none of the above constraints holds. The bold edges show the paths in the tree corresponding to the data point $x=(0,0)$. In this case, $x$ falls into a single bucket $w_{1}$.

Theorem 2

No decision tree exists with overall complexity below $O(N^{\lambda^{*}})$.

Theorem 3

The decision tree construction of Algorithm 4 described in Remark 7 results in a tree with complexity $O(N^{\lambda^{*}})$.

In other words, Theorem 3 proves that the tree $G$ constructed by Algorithm 4 satisfies ((34)-(38)), and Theorem 2 shows that this is the optimal decision tree. For proofs of Theorems 2 and 3, see Sections 6.1 and 6.2. Note that, not only Theorem 3 guarantees that the number of the nodes in our decision tree is bounded by $O(N^{\lambda^{*}})$ but also it guarantees that the run time for mapping the data points to the decision tree and the number of comparisons that we need to do for the nodes with the collision is bounded by $O(N^{\lambda^{*}})$, see complexity equation (33) in Section 4.1.

Theorem 4** (Noise Robustness)**

Assume that the decision tree of Algorithm 4 described in (57) is constructed based on distribution $P=[p_{ij}]$ while the data is generated from unknown $P^{\prime}=[p^{\prime}_{ij}]$ satisfying $\frac{p_{ij}}{1+\epsilon}\leq p^{\prime}_{ij}\leq{p_{ij}}{(1+\epsilon)}$ for all $1\leq i\leq k,1\leq j\leq l$ and some $\epsilon>0$. Then, it is possible to design an algorithm with complexity $O(N^{\lambda^{*}(p)+3c_{d}\log(1+\epsilon)})$ where $c_{d}=\frac{{(\lambda^{*}-\min(1,\delta))}}{\log(\max_{i,j}\min(\frac{p_{ij}}{p_{i}^{\mathcal{A}}},\frac{p_{ij}}{p_{j}^{\mathcal{B}}}))}$ while maintaining arbitrary high true positive rate $TP$.

In other words, Theorem 4 proves that the tree $G$ constructed by Algorithm 4 is noise robust, i.e., for the case where our understanding from the distribution is not exact, the arbitrary high true positive rate $TP$ can still be maintained while the complexity increases linearly with the noise. For proof of Theorem 4, see Section 6.3. For a high level intuition behind the relationship between Algorithm 4 and Theorems 2 and 3, note that: Every constraint in the decision tree defined by constraints (57) directly relates to the constraints $(\eqref{order1}-\eqref{order3})$. Therefore, we expect this decision tree to be optimal.

Here, we present an intuition behind proof of Theorem 2. For rigorous proof of Theorem 2, see Section C.

6.1 Intuition behind proof of Theorem 2

Note that, from definition of $\Phi(v)$, ${\Psi}^{\mathcal{A}}(v)$, ${\Psi}^{\mathcal{B}}(v)$ and $\Psi(v)$ in $(\eqref{def_a}-\eqref{def_b})$, we have $\Phi(w_{ij})=p_{ij}\Phi(v)$ where $v$ is the parent of $w_{ij}$. Similar equations hold for ${\Psi}^{\mathcal{A}}(v)$, ${\Psi}^{\mathcal{B}}(v)$ and $\Psi(v)$. On the other hand, from (52), note that for any $(\mu,\nu,\eta)\in\mathcal{I}$ we have $\sum_{1\leq i\leq k,1\leq j\leq l}p_{ij}^{1+\mu+\nu-\eta}{(p_{i}^{\mathcal{A}})}^{-\mu}{(p_{j}^{\mathcal{B}})}^{-\nu}=1$. Therefore, we expect to have similar relation between $\Phi(v)$, ${\Psi}^{\mathcal{A}}(v)$, ${\Psi}^{\mathcal{B}}(v)$, $\Psi(v)$, i.e.,

[TABLE]

for any $(\mu,\nu,\eta)\in\mathcal{I}$. On the other hand, from $(\eqref{alphaav}-\eqref{gammydef})$, we have $\alpha(G)=\sum_{v\in V_{Buckets}(G)}\Phi(v)$. Similar definitions hold for $\beta(G)$, $\gamma^{\mathcal{A}}(G)$ and $\gamma^{\mathcal{B}}(G)$. Therefore, from convexity of the function $f(\theta,\theta_{1},\theta_{2},\theta_{3})={\theta}^{1+\rho_{1}+\rho_{2}+\rho_{3}}{\theta_{1}}^{-\rho_{1}}{\theta_{2}}^{-\rho_{2}}{\theta_{3}}^{-\rho_{3}}$ we conclude

[TABLE]

using the lower bounds on $\frac{\alpha(G)}{\beta(G)}$, $\frac{\alpha(G)}{\gamma^{\mathcal{A}}(G)}$, $\frac{\alpha(G)}{\gamma^{\mathcal{B}}(G)}$ and ${\alpha(G)}$ in $(\eqref{order1}-\eqref{order3})$ and (59), we conclude $\lambda\geq\lambda^{*}$.

Here, we present an intuition behind proof of Theorem 3. For rigorous proof of Theorem 3, see Appendix D.

6.2 Intuition behind proof of Theorem 3

Here, our goal is to prove that the tree construction steps in Algorithm 4, i.e., (57) result in a tree that satisfies

There is at least one node in the tree which is accepted as a bucket. 2. 2.

$(\eqref{order111}-\eqref{order3})$ holds. 3. 3.

Expected number of nodes in the tree is bounded by $O(N^{\lambda^{*}})$.

Let us intuitively prove all these three statements one by one.

Consider a node with the depth equal to $n^{*}$ given in (55). Moreover, assume that the number of times we have $a_{i}$ at $s$-th position of $x$ and $b_{j}$ at $s$-th position of $y$ for all $1\leq s\leq n^{*}$ are $n^{*}r^{*}_{ij}$. Then, we argue that this node is not pruned and accepted as a bucket. In order to understand why intuitively it is true, we verify that this node is accepted as a bucket, i.e., we have to verify that $\frac{\Phi(v)}{\Psi(v)}\geq{N^{1+\delta-\lambda^{*}}}{p_{0}q_{0}}$. This is true as

[TABLE]

(60) follows from (21), (22) and the definition of node $v$, i.e., the number of times we have $a_{i}$ at $s$-th position of $x$ and $b_{j}$ at $s$-th position of $y$ are $n^{*}r^{*}_{ij}$. (61) is concluded from the definition of $r^{*}_{ij}$ in (54). For further details, see Section D. 2. 2.

Let us prove (35) as the rest of $(\eqref{order112}-\eqref{order3})$ follow similarly. From Algorithm 4, for all the buckets we have

[TABLE]

Note that, at least there is one node that is accepted a a bucket. On the other hand, $\alpha(G)=\sum_{v\in V_{Buckets}(G)}\Phi(v)$ and $\beta(G)=\sum_{v\in V_{Buckets}(G)}\Psi(v)$. Therefore, we conclude that

[TABLE]

Note that $\frac{\sum_{i}{a_{i}}}{\sum_{i}b_{i}}\geq c$ if $\frac{a_{i}}{b_{i}}\geq c$ and $b_{i}>0$ for any $i$. 3. 3.

Finally, we prove that number of nodes in the tree is bounded by $O(N^{\lambda^{*}})$. Note that $\sum_{i,j}p_{ij}=1$, therefore, we expect to have $\sum_{v\in V_{l}(G)}\Phi(v)=1$. On the other hand, $\Phi(v)$ is expected to be greater than $N^{-\lambda^{*}}$ for the intermediate nodes from (57). Therefore, it is concluded that $|V_{l}(G)|$ is at most $N^{\lambda^{*}}$. This results in $|V(G)|=O(N^{\lambda^{*}})$ as for any decision tree we have $|V(G)|\leq 2|V_{l}(G)|$. For details, see Appendix D.

.

Here, we present an intuition behind proof of Theorem 4. For rigorous proof of Theorem 4, see Appendix E.

6.3 Intuition behind proof of Theorem 4

Define $\Phi^{\prime}(v),\alpha^{\prime}(G),\beta^{\prime}(G),\#bands^{\prime}$ for $p^{\prime}_{ij}$ the same way as $\Phi(v),\alpha(G),\beta(G),\#bands$ for $p_{ij}$. As $\frac{p_{ij}}{1+\epsilon}\leq p^{\prime}_{ij}\leq{p_{ij}}{(1+\epsilon)}$, we expect $\Phi(v)$ to be bounded as

[TABLE]

Therefore, from $\eqref{alphaav}$ we have

[TABLE]

We conclude similar inequalities for $\gamma^{\mathcal{A}}(G)$ and $\gamma^{\mathcal{B}}(G)$, while for $\beta(G)$ we have

[TABLE]

From (51), $\#bands$ is inversely related to $\alpha(G)$. Therefore, $\#bands^{\prime}$ can be bounded from above by ${(1+\epsilon)}^{d}\#bands$. As a result, from (33) the total complexity is bounded by ${(1+\epsilon)}^{3d}N^{\lambda^{*}}$ from above. Assume that $d$ is bounded by $c_{d}\log N$ for some constant $c_{d}$. Therefore, we conclude Theorem 4. In order to see why $d\leq c_{d}\log N$, note that for all the leaf nodes, we have $\frac{\Phi(v)}{{\Psi}^{\mathcal{A}}(v)}\leq{N^{1+\delta-\lambda^{*}}}\max_{i,j}\frac{p_{ij}}{p_{i}^{\mathcal{A}}}p_{0}q_{0}$. On the other hand, $\frac{\Phi(v)}{{\Psi}^{\mathcal{A}}(v)}$ can be bounded by {\big{(}\min_{i,j}\frac{p_{ij}}{p_{i}^{\mathcal{A}}}\big{)}}^{d} from below. Therefore, we conclude that $d=c_{d}\log N$ for some constant $c_{d}$.

7 Experiments and Observations

Experiment 1

In this experiment, we compare the complexity of ForestDSH with the algorithm proposed by Dubiner in [41]. Here, we set $\mathcal{A}=\mathcal{B}=\{0,1\},S=1000,M=N=1000,000$. For ForestDSH, We use Algorithm 5 to derive $\mu^{*},\nu^{*},\eta^{*},\lambda^{*}$ and $\delta$ while training $C_{1}$, $C_{2}$ and $C_{3}$ to get the best complexity in Algorithm 4. For Dubiner’s algorithm, equation $(126)$ in [41] is computationally intractable which makes it impossible to compute the complexity for the algorithm presented there for general probability distributions. However, in the special case where $P(p)=\begin{bmatrix}\frac{p}{2}&\frac{1-p}{2}\\ \frac{1-p}{2}&\frac{p}{2}\end{bmatrix}$, $0.5\leq p\leq 1$ (hamming distance), a solution has been provided for computing the complexity in [41], and we implemented that solution (see Appendix F for the detail of implementation), and compared it to ForestDSH. Note that currently no source code for the implementation of Dubiner algorithm is available. Figure 5, shows that Dubiner algorithm’s performance is worse than that of ForestDSH.

Observation 1

In this observation, we reformulate the problem of solving (2) to the minimum inner product search problem (MIPS) [35] by transferring the data points from $\mathcal{A}^{S}$ and $\mathcal{B}^{S}$ to $R^{klS}$ in a way that \log\big{(}\frac{\mathbb{P}(x,y)}{\mathbb{Q}(x,y)}\big{)} is equal to the dot product in this new space. We transformed $x\in\mathcal{A}^{S}$ and $y\in\mathcal{B}^{S}$ to $T(x)\in\mathbb{R}^{klS}$ and $T(y)\in\mathbb{R}^{klS}$ as follows:

[TABLE]

Then, we have $\log\left(\frac{\mathbb{P}(x,y)}{\mathbb{Q}(x,y)}\right)=<T(x),T(y)>$ where $<.,.>$ stands for the inner product in $\mathbb{R}^{klS}$. In other words, given any $x=(x_{1},\cdots,x_{S})$, for each $x_{s}$ we transform it to a $kl\times 1$ vector with $l$ non-zero elements. Similarly, given any $y=(y_{1},\cdots,y_{S})$, each $y_{s}$ is transformed to a $kl\times 1$ vector with $k$ non-zero elements. Therefore, finding pairs of data points with large $\frac{\mathbb{P}(x,y)}{\mathbb{Q}(x,y)}$ is equivalent to finding transformed data points with large dot product. Using this transformation, in Appendix G we show that the angle between both the true pairs and false pairs will be nearly $\frac{\pi}{2}$ for almost all the probability distributions ($\frac{S_{0}}{M^{2}}\approx 0$, using the notation from [35]). It is well known that MIPS performs poorly in detection of pairs that are nearly orthogonal [35]. Therefore, solving (2) by transforming it to a MIPS problem and using existing approaches fails. Note that, [35] is based on data independent hashes for Euclidean distance that are not the state of the art. Recently, data dependent hashes have been introduced for Euclidean distance that improve on their data independent counterparts [33, 34], while currently there is no data dependent strategy for maximum inner product search. Therefore, MIPS is currently unable to solve (1) efficiently.

Experiment 2

In this experiment, we compared the complexity for the three algorithms LSH-hamming, MinHash and ForestDSH for a range of probability distributions. We benchmark the three methods using matrices ${P}(t)={P}_{1}(1-t)+{P}_{2}t$ where $0\leq t\leq 1$, ${P}_{1}=\begin{bmatrix}0.345&0\\ 0.31&0.345\end{bmatrix}$, ${P}_{2}=\begin{bmatrix}0.019625&0\\ 0.036875&0.9435\end{bmatrix}$, and $\delta=1$, i.e., $M=N$. The selection of ${P}_{1}$ was such that the complexity of MinHash minus the complexity of LSH-hamming was maximized. ${P}_{2}$ was selected such that the complexity of LSH-hamming minus the complexity of MinHash was maximized. Fig. 6 $(a)$ shows the theoretical complexities of MinHash, LSH-hamming and ForestDSH for each matrix. See Appendix H, for the details on the derivation of complexities for MinHash, LSH-hamming and ForestDSH. For instance, for ${P}_{1}=\begin{bmatrix}0.345&0\\ 0.31&0.345\end{bmatrix}$, the theoretical per query complexities of MinHash, LSH-hamming and ForestDSH are equal to $0.5207$, $0.4672$ and $0.4384$, respectively. We further consider $N$ data points of dimension $S$, $\{x^{1},\cdots,x^{N}\}$ and $\{y^{1},\cdots,y^{N}\}$ where each $(x^{i},y^{i})$ is generated from ${P}(t)$, and $x^{i}$ is independent from $y^{j}$ for $i\neq j$ $(N=2000,S=2000)$. Then, we used ForestDSH, LSH-hamming and MinHash to find the matched pairs. In each case, we tuned $\#rows$666In MinHash and LSH-Hamming, we start with $\#bands\times\#rows$ randomly selected hashes from the family of distribution sensitive hashes, where $\#rows$ is the number of rows and $\#bands$ is the number of bands. We recall a pair $(x,y)$, if $x$ and $y$ are hashed to the same value in all $\#rows$ rows and in at least one of the $\#bands$ bands. and $\#bands$ to achieve $99\%$ true positive (recall) rate. Total simulation time for each of the three methods is plotted for each probability distribution in Figure 6 $(b)$. The simulation times in Figure 6 $(b)$ are consistent with the theoretical guarantees in Figure 6 $(a)$. Figure 6 parts $(b)$ and $(c)$ show that for sparse matrices, $(t\approx 1)$, MinHash and ForestDSH outperform LSH-hamming. In denser cases, $(t\approx 0)$, LSH-hamming and ForestDSH outperform MinHash. For $(t\leq 0.4)$, ForestDSH outperforms both MinHash and LSH-hamming. Note that, for the case when the data is coming from sparse distribution, MinHash beats ForestDSH in practice (as opposed to Theory). This is because for sparse data the total complexity tends to its minimum, i.e., $O(N)$. ForestDSH is inefficient compared to MinHash for small number of $N$ and when the data is sparse as the constant terms and $logN$ terms which play a significant role in this case are not optimized in ForestDSH. In Figure 7, we further plotted $V(G(N))$, $\frac{\alpha(G(N))}{\beta(G(N))}$, $\frac{\alpha(G(N))}{\gamma^{\mathcal{A}}(G(N))}$ and $\frac{\alpha(G(N))}{\gamma^{\mathcal{B}}(G(N))}$ as a function of $N$ for trees $G(N)$ constructed by Algorithm 4 for ${P}(t=0.25)$, where $M=N$. As predicted by Theorem 3, we observed that these quantities grow/decay proportional to $N^{\lambda^{*}}$, $N^{1+\delta-\lambda^{*}}$, $N^{1-\lambda^{*}}$ and $N^{\delta-\lambda^{*}}$, respectively. In Figure 8, true positive rate and total complexity are plotted in terms of $\#bands$ for the case of $N=20,000$, $S=10,000$ and the probability distribution ${P}_{2}$. From (50), i.e., $TP=1-{(1-\alpha(G))}^{\#bands}$, we expect true positive rate to increase while $\#bands$ increases (see Figure 8 $(a)$). On the other hand, from (51) and (33), total complexity $(N^{\lambda^{*}})$ is expected to linearly increase while $\#bands$ increases (see Figure 8 $(b)$).

Experiment 3

In this experiment, we examine the robustness of our algorithm from Theorem 4 for matrices considered in Experiment 2 while $M=N=300,000$. For each ${P}(t)$, we generated random matrices ${P}^{\prime}(t)$ within $1+\epsilon$ of ${P}(t)$, see Theorem 4. We derive the practical complexity and theoretical complexity, i.e., $\lambda^{*}$ from Definition 4. Moreover, for each of these matrices, we derive the worst case complexity in the case of $\epsilon=0.03$. This is sketched in Figure 9.

Experiment 4

In this experiment, we used the mass spectral data from human microbiome [43], extracted using MSCluster [1]. In this problem, there are $N$ spectra $X={x^{1},\cdots,x^{N}}$, and a query spectrum $y$. Our goal is to find the spectra $x^{i}$ that maximize a probabilistic model $\mathbb{P}(y|x^{i})$. $\mathbb{P}(y|x)$ is learned assuming that it can be factorized to i.i.d. components. To learn $p(y\mid x)$, each mass spectra is sorted based on the peak intensities, and in order to reduce the number of parameters we need to learn, instead of the peak ranks we use log of the peak ranks 777$\log_{b}Rank$ of a peak is defined as the $\log_{b}$ of its $rank$. For any natural number $n$, $\log_{b}Rank=n$ for the peaks at rank $\{b^{n-1},\cdots,b^{n}-1\}$, e.g., $\log_{2}Rank(m)=3$ for $m\in\{4,5,6,7\}$. Joint probability distribution of logRanks for the data from [1] is shown in Figures 10 $(a)$ ($4\times 4$ data matrix is obtained using $\log_{4}Rank$), $(b)$ ($8\times 8$ data matrix is obtained using $\log_{2}Rank$) and $(c)$ ($51\times 51$ data matrix is obtained not using any $\log_{b}Rank$).. Using these data, we learn the joint probability distribution $p(\log Rank(y_{s})=i\mid\log Rank(x_{s})=j)$. Among 90753 data points, we used 70753 for training, and 20000 for the test. The joint probability distribution of matching pairs are shown in Figures 10, before and after logRank transformation.

*After learning this probabilistic model, we applied ForestDSH method to the test mass spectra (after $\log Rank$ transformation), and we were able to speed up the search of $20000$ mass spectra nine times, in comparison to brute force search while maintaining true positive rate of $90\%$. For ForestDSH method, the total runtime is $705,197ms$, while for brute force, the total runtime is 6,022,373ms (eight times slower than ForestDSH), and for LSH, the runtime for achieving $90\%$ TP rate is 5,495,518ms (seven times slower than ForestDSH). The amount of memory used peaks at 220MB. *

Experiment 5

In this experiment, we considered a set of $50,000$ images ($25,000$ pairs of images) from CIFAR-10 database [44]. Each pair of images consists of a grey-scale $32*32$ pixels image and a noisy version of the image constructed by adding independent Gaussian noise $\mathcal{N}(\mu=0.5,\sigma=0.039)$ and discretizing pixels to binary. ForestDSH is able to detect true pairs with the success rate of $99\%$ in 1,970 seconds while brute force approach detects true pairs in 19,123 seconds (nine times slower than ForestDSH) and LSH detects true pairs within 15,080 seconds (seven times slower than ForestDSH).

7.1 Codes

For the codes, see https://github.com/mohimanilab/ForestDSH.

8 Conclusion

ForestDSH algorithm proposed in this paper is comprehensive and efficient in the sense that for a wide range of probability distributions it enables us to capture the difference between pairs coming from a joint probability distributions and independent pairs. This algorithm is built upon a family of distribution sensitive hashes and is designed using a decision tree structure which is constructed recursively. Moreover, we prove that the decision tree introduced here has a complexity of $O(N^{\lambda^{*}})$ and there is no decision tree with overall complexity below it. We prove that this algorithm outperforms existing state of the art approaches in specific range of distributions, in theory and practice and enabled speeding up the spectral library search in mass spectrometry by a factor of nine. Finally, we should note that ForestDSH has some limitations which are discussed as follows. ForestDSH assumes that the probability distribution function is factorizable to i.i.d. components. Generalizing ForestDSH to the case of general probability distribution function is an open problem and our results here open a path towards solving the general probability mass function problem. Moreover, the special case of factorizable models are simple but crucial models that have been widely used in computational biology, e.g., see [3], and other areas of data sciences. In the future, we will address the more general problem of Markov chain models. Note that ForestDSH performs faster than MinHash if the data is not sparse. While for a small number of datapoints MinHash is superior to ForestDSH on sparse data, for a larger number of datapoints ForestDSH is performing the same as MinHash. This is mainly due to the fact that constant and logarithmic terms vanish in ratio as the number of datapoints grows.

Appendix A Proof of Lemmas 1 and 2

A.1 Proof of Lemma 1.

Lemma 1 is proved by induction on the depth of the node $v$. For example, consider $v$ and its children $w_{ij}=f(v,a_{i},b_{j})$. From (21), we have $\Phi(w_{ij})=\Phi(v)p_{ij}$. Therefore, (25) is proved by induction as follows. Assume (25) holds for any node with depth less than $d$. Consider the node $w_{ij}$ which is a children of $v$, i.e., $w_{ij}=f(v,a_{i},b_{j})$ and $depth(w_{ij})=d$.

[TABLE]

See Section 4 for the Definition of $perm_{z}$. Note that ${({perm_{z}}(x))}_{d}$ stands for $d$-th entry of the vector ${({perm_{z}}(x))}$. (A.1) follows from induction assumption for the nodes with depth less than $d$, (LABEL:tp3) is a result of i.i.d. assumption (2), i.e., $\mathbb{P}(y\mid x)=\prod_{s=1}^{S}p(y_{s}\mid x_{s})$ and the defnition of $H_{z}^{\mathcal{A}}(x)$ in (19). $(\eqref{famm1}-\eqref{fam2})$ follow similarly.

A.2 Proof of Lemma 2.

Using $(\eqref{fam1}-\eqref{fam2})$, constraints $(\eqref{family1}-\eqref{family2})$ hold for $\alpha=\alpha(G)$, $\beta=\beta(G)$, $\gamma^{\mathcal{A}}=\gamma^{\mathcal{A}}(G)$ and $\gamma^{\mathcal{B}}=\gamma^{\mathcal{B}}(G)$. Since no two buckets are ancestor/descendant of each other, we have

[TABLE]

Therefore, (15) holds. This completes the proof that $H_{z}^{\mathcal{A}}(x)$ and $H_{z}^{\mathcal{B}}(y)$ defined in (19) and (20) are

$(\alpha(G),\beta(G),\gamma^{\mathcal{A}}(G),\gamma^{\mathcal{B}}(G))$-sensitive.

Appendix B Deriving $\mu^{*}$, $\nu^{*}$, $\eta^{*}$, $\lambda^{*}$, $p_{0}$, $q_{0}$ and $\delta$ for $\mathbb{P}$, $M$ and $N$

In this section, an algorithm for deriving $\mu^{*}$, $\nu^{*}$, $\eta^{*}$, $\lambda^{*}$, $p_{0}$, $q_{0}$ and $\delta$ for $\mathbb{P}$, $M$ and $N$ for the probability distribution $\mathbb{P}$ and $\delta$ is presented for a given probability distribution $\mathbb{P}$.

Remark 8

In Algorithm 5, the parameters $\mu^{*}$, $\nu^{*}$, $\eta^{*}$ and $\lambda^{*}$ could be derived from newton method too.

Appendix C Proof of Theorem 2

In order to prove Theorem 2, we first state the following two lemmas.

Lemma 3

The function $f(\theta,\theta_{1},\theta_{2},\theta_{3})={\theta}^{1+\rho_{1}+\rho_{2}+\rho_{3}}{\theta_{1}}^{-\rho_{1}}{\theta_{2}}^{-\rho_{2}}{\theta_{3}}^{-\rho_{3}}$ is a convex function on the region $(\theta,\theta_{1},\theta_{2},\theta_{3})\in\mathbb{R}^{4+}$ where $(\rho_{1},\rho_{2},\rho_{3})\in\mathbb{R}^{3+}$. 888For any natural number $n$, $\mathbb{R}^{n+}$ denotes as the set of all $n$-tuples non-negative real numbers.

Lemma 4

\sum_{v\in V_{Buckets}(G)}{\big{(}\Phi(v)\big{)}}^{1+\mu+\nu-\eta}{\big{(}{\Psi}^{\mathcal{A}}(v)\big{)}}^{-\mu+\eta}{\big{(}{\Psi}^{\mathcal{B}}(v)\big{)}}^{-\nu+\eta}{\big{(}\Psi(v)\big{)}}^{-\eta}\leq 1* for any $(\mu,\nu,\eta)\in\mathcal{I}$.*

Proof of Lemma 3 and Lemma 4, are relegated to Appendices C.1 and C.2, respectively. Consider $(\mu^{*},\nu^{*},\eta^{*})$ that satisfy (53). For any decision tree satisfying $(\eqref{order1}-\eqref{order3})$, we have:

[TABLE]

where (79) holds to the convexity of $f(\theta,\theta_{1},\theta_{2},\theta_{3})={\theta}^{1+\rho_{1}+\rho_{2}+\rho_{3}}{\theta_{1}}^{-\rho_{1}}{\theta_{2}}^{-\rho_{2}}{\theta_{3}}^{-\rho_{3}}$ in Lemma 3 and (80) follows from Lemma 4. Therefore, we have

[TABLE]

On the other hand, using (80) and the definitions of $\alpha(G)$, $\beta(G)$, $\gamma^{\mathcal{A}}(G)$ and $\gamma^{\mathcal{B}}(G)$ in $(\eqref{alphdef}-\eqref{gammydef})$, we have

[TABLE]

Therefore, from $(\eqref{order1}-\eqref{order3})$ and (82) we have

[TABLE]

Therefore, we conclude Theorem 2 as follows

[TABLE]

C.1 Proof of Lemma 3

The Hessian matrix for $f(\theta,\theta_{1},\theta_{2},\theta_{3})$ is represented as

[TABLE]

where $W=\frac{\rho_{1}+\rho_{2}+\rho_{3}}{\theta}$, $W_{i}=\frac{\rho_{i}}{\theta_{i}}$ for any $i\in\{1,2,3\}$. In order to show that the function $f(\theta,\theta_{1},\theta_{2},\theta_{3})$ is a convex function it is necessary and sufficient to prove that $H(\theta,\theta_{1},\theta_{2},\theta_{3})$ is positive semidefinite on $\mathbb{R}^{4+}$ On the other hand, for positive semidefinite matrices we have

For any non-negative scalar $a$ and positive semidefinite matrix $M$, $aM$ is positive semidefinite. 2. 2.

For positive semidefinite matrices $M_{1}$ and $M_{2}$, $M_{1}+M_{2}$ is positive semidefinite.

As $f(\theta,\theta_{1},\theta_{2},\theta_{3})>0$ for any $\theta,\theta_{1},\theta_{2},\theta_{3}$, it is sufficient to prove that $\frac{H(\theta,\theta_{1},\theta_{2},\theta_{3})}{f(\theta,\theta_{1},\theta_{2},\theta_{3})}$ is positive semidefinite. Define, $M_{1}=\begin{bmatrix}W^{2}&WW_{1}&WW_{2}&WW_{3}\\ W_{1}W&W_{1}^{2}&W_{1}W_{2}&W_{1}W_{3}\\ W_{2}W&W_{2}W_{1}&W_{2}^{2}&W_{2}W_{3}\\ W_{3}W&W_{3}W_{1}&W_{3}W_{2}&W_{3}^{2}\end{bmatrix}$ and $M_{2}=\begin{bmatrix}\frac{\rho_{1}+\rho_{2}+\rho_{3}}{\theta^{2}}&-\frac{\rho_{1}}{\theta\theta_{1}}&-\frac{\rho_{2}}{\theta\theta_{2}}&-\frac{\rho_{3}}{\theta\theta_{3}}\\ -\frac{\rho_{1}}{\theta\theta_{1}}&\frac{\rho_{1}}{\theta_{1}^{2}}&0&0\\ -\frac{\rho_{2}}{\theta\theta_{2}}&0&\frac{\rho_{2}}{\theta_{2}^{2}}&0\\ -\frac{\rho_{3}}{\theta\theta_{3}}&0&0&\frac{\rho_{3}}{\theta_{3}^{2}}\end{bmatrix}$. Therefore, we have $M_{1}+M_{2}=\frac{H(\theta,\theta_{1},\theta_{2},\theta_{3})}{f(\theta,\theta_{1},\theta_{2},\theta_{3})}$. In order to prove that $f(\theta,\theta_{1},\theta_{2},\theta_{3})$ is positive semidefinite, it is sufficient to prove that $M_{1}$ and $M_{2}$ are positive semidefinites. The matrices $M_{1}$ and $M_{2}$ are positive semidefinite as for any non-zero vector $z=\begin{bmatrix}a&b&c&d\end{bmatrix}$, we have $zM_{1}z^{T}\geq 0$ and $zM_{2}z^{T}\geq 0$, i.e.,

[TABLE]

where (91) is concluded as $\rho_{1},\rho_{2},\rho_{3}\geq 0$.

C.2 Proof of Lemma 4

First of all, note that $\Psi(v)={\Psi}^{\mathcal{A}}(v){\Psi}^{\mathcal{B}}(v)$. Let us define

[TABLE]

We show

[TABLE]

by induction on the number of nodes in the tree. If the tree has only one node, i.e., root, then (93) is holds as $\Phi({root})=1$, ${\Psi}^{\mathcal{A}}({root})=1$ and ${\Psi}^{\mathcal{B}}({root})=1$ from the definition of $\Phi(v)$, ${\Psi}^{\mathcal{A}}(v)$ and ${\Psi}^{\mathcal{B}}(v)$ in $(\eqref{def_a}-\eqref{def_a3})$. Assume that (93) holds for any decision tree with $|G|<Z$. Our goal is to prove that (93) holds for a decision tree with $|G|=Z$. Assume $v_{11}$ is the node with maximum length in $G$ and consider a tree $G^{\prime}$ constructed by removing $v_{11}$ and all its siblings $v_{ij},1\leq i\leq k,1\leq j\leq l$ belonging to the same parent $v$. In other words, for the tree $G^{\prime}$ we have999Note that, $V_{b}(G^{\prime})$ satisfies bucket-list property, e.g., there is no bucket in the tree that is ancestor of another bucket.

[TABLE]

(95) is true as for the gragh $G^{\prime}$, the node $v$ in now a leaf node while the nodes $v_{11},\cdots,v_{kl}$ are removed. Then, we have

[TABLE]

where (98) holds from the definition of tree $G^{\prime}$, (98) follows from the recursive definition of $\Phi(v)$, ${\Psi}^{\mathcal{A}}(v)$ and ${\Psi}^{\mathcal{B}}(v)$ in $(\eqref{def_a}-\eqref{def_a3})$ and (99) holds, note that from the definition of $\mu$, $\nu$ and $\eta$ in (52), i.e.,

[TABLE]

Therefore, we conclude that

[TABLE]

Note that, the inequality (98) becomes an equality only in cases where the tree is homogeneous and none of the children are pruned.

Appendix D Proof of Theorem 3

In order to prove Theorem 3, we first present the following lemma.

Lemma 5

Given a fixed $N\in\mathbb{N}$ and probability distribution $\mathbb{P}$, consider the following region $R_{\lambda}$

[TABLE]

101010Recall that, for simplicity we use the notation $\sum_{i,j}$ and $\prod_{i,j}$ instead of $\sum_{1\leq i\leq k,1\leq j\leq l}$ and $\prod_{1\leq i\leq k,1\leq j\leq l}$, respectively.

Then, $(\lambda^{*},r_{ij}^{*},n^{*})$ defined in Definition 4 is a member of $\mathcal{R}(\lambda,r_{ij},n)$.

The proof of Lemma 5 is relegated to Appendix D.4. Let us prove that the following tree construction steps in Algorithm 4 result in a tree that satisfies $(\eqref{order1}-\eqref{order3})$.

[TABLE]

Consider the set of $r^{*}_{ij}=p_{ij}^{1+\mu^{*}+\nu^{*}-\eta^{*}}{(p_{i}^{\mathcal{A}})}^{-\mu^{*}}{(p_{j}^{\mathcal{B}})}^{-\nu^{*}}$ and $n^{*}=\frac{(\max(1,\delta)-\lambda^{*})\log N}{\sum{r_{ij}^{*}\log\frac{p_{ij}}{r_{ij}^{*}}}}$. Note that we assume $p_{ij}$ and $q_{ij}$ are non-zero111111Note that in the cases where $q_{ij}$ is zero, then from the definition of $q_{ij}$, $p_{ij}$ would also be equal to zero. Therefore, we will ignore those branches during the tree construction.. Consider $n_{ij}=\lceil n^{*}r_{ij}^{*}\rceil$ if $r^{*}_{ij}>\frac{1}{2}$ and $n_{ij}=\lfloor n^{*}r^{*}_{ij}\rfloor$ if $r^{*}_{ij}\leq\frac{1}{2}$. Therefore, we have $n^{*}-kl<\sum_{ij}n_{ij}\leq n^{*}$. For any $v\in V(G)$, define the set ${S}_{ij}(v)$ as follows

[TABLE]

where $depth(v)$ is the depth of node $v$ in the tree, $Seq^{\mathcal{A}}_{s}(v)$ and $Seq^{\mathcal{B}}_{s}(v)$ stand for position $s$ in the strings $Seq^{\mathcal{A}}(v)$ and $Seq^{\mathcal{B}}$(v), respectively. Now, consider a node $v^{*}$ in the graph that satisfies the following constraints:

[TABLE]

The number of nodes $v$ that satisfy this constraint is $\binom{n^{*}}{n_{11},\cdots,n_{kl}}$. Moreover, define

[TABLE]

D.1 Node $v^{*}$ or one of its ancestors is designated as a bucket by Algorithm 4

Here, we prove that the node $v$ or one of its ancestors is designated as a bucket by Algorithm 4. In order to show this, we need to prove that:

[TABLE]

where $p_{0}$ and $q_{0}$ are defined as $\prod_{i,j}p_{ij}$ and $\min(\prod_{i,j}q_{ij},\prod_{i}{(p_{i}^{\mathcal{A}})}^{l},\prod_{j}{(p_{j}^{\mathcal{B}})}^{k})$. Note that, $\Phi(v^{*})$, $\Psi(v^{*})$, ${\Psi}^{\mathcal{A}}(v^{*})$ and ${\Psi}^{\mathcal{B}}(v^{*})$ are computed as follows

[TABLE]

Therefore, from Lemma 5 and $(\eqref{..1}-\eqref{..4})$ we conclude $(\eqref{qmmm1}-\eqref{qttttt2})$. This means $v^{*}$ or one of its ancestors is an accepted bucket.

D.2 Proof of bounds ((35)-(38))

First of all, we derive a lower bound on $\alpha(G)$ as follows.

[TABLE]

where $V_{n_{11},\cdots,n_{kl}}$ is the set of nodes that satisfies (110). $|V_{n_{11},\cdots,n_{kl}}|$ is lower bounded as

[TABLE]

for some constant $c=\frac{(2\pi)^{\frac{1-kl}{2}}e^{-kl}}{(kl)!}{(n^{*})}^{\frac{1-kl}{2}}$ depending on $n$, $k$ and $l$. (122) is true as for any natural number $m$ we have $\sqrt{2\pi m}{(\frac{m}{e})}^{m}\leq m!<\sqrt{2\pi m}{(\frac{m}{e})}^{m}e$. (125) follows as $a\log\frac{1}{a}$ is an increasing function for $0\leq x\leq 0.5$, and a decreasing function for $0.5\leq x\leq 1$. Therefore, from (121) and (125), ${\alpha(G)}\geq N^{\max(1,\delta)-\lambda^{*}}$ is concluded. Similarly, $(\eqref{order1}-\eqref{order3})$ are proved as follows.

[TABLE]

where (126) and (128) is concluded from $(\eqref{i-i=1}-\eqref{iii2})$ and the fact that $\frac{\sum_{i}{a_{i}}}{\sum_{i}b_{i}}\geq c$ if $\frac{a_{i}}{b_{i}}\geq c$ and $b_{i}>0$ for any $i$.

D.3 Bounding number of nodes in the tree, i.e., (34)

The number of nodes in the decision tree defined in (57), is bounded by $O(N^{\lambda^{*}})$ using the following three lemmas.

Lemma 6

For any leaf node $v$ of the decision tree defined in (57), we have

[TABLE]

Lemma 7

For any tree $G$, the summation of $\Phi(v)$ over all the leaf nodes is equal to one, i.e., $\sum_{v\in V_{l}}\Phi(v)=1$.

Lemma 8

The number of nodes in the decision tree defined in (57) is at most two times of the number of leaf nodes.

For proof of Lemmas 6, 7 and 8, see Appendices D.5, D.6 and D.7. Therefore, we have

[TABLE]

where (130) follows from Lemma 7, (131) is true from (167) and (133) is concluded from Lemma 8. Therefore, we conclude that $|V(G)|=O(N^{\lambda^{*}})$.

D.4 Proof of Lemma 5

Consider the optimization problem of finding the member of $(\lambda,r_{ij},n)\in\mathcal{R}(\lambda,r_{ij},n)$ with minimum $\lambda$. This optimization problem is a convex optimization problem121212Note that $f_{1}(n)=\frac{1}{n}$, $f_{2}(r)=r\log r$ and $f_{3}(r)=ar$ are convex functions. Therefore, as the sum of the convex functions is a convex function, the optimization problem $(\eqref{l0}-\eqref{l---1})$ is in the form of convex optimization problem

$\displaystyle\underset{\lambda}{\mbox{Minimize~{}}}~{}f(x),$

(134)

$\displaystyle{\mbox{subject to~{}}}~{}\{g_{i}(x)$ $\displaystyle\leq$ $\displaystyle 0,i\in\{1,\cdots,m\},$

(135)

$\displaystyle h_{j}(x)$ $\displaystyle=$ $\displaystyle 0,j\in\{1,\cdots,p\}\},$

(136)

where $x\in\mathbb{R}^{n}$, $f(x)$ and $g_{i}(x)$ are convex functions and $h_{j}(x)$ is affine functions.. Therefore, writing the KKT conditions, we have

[TABLE]

where

[TABLE]

From (138), $\mu_{1ij}$ is zero if $r_{ij}$ is a non-zero number. Therefore, we only keep $i$ and $j$ where $r_{ij}\neq 0$ and $\mu_{1ij}=0$.

[TABLE]

Consider the following two cases.

$\mu_{2}=0$. In this case, all the constraints are affine functions and therefore we have a linear programming problem and the feasible set of this linear programming problem is a polyhedron. From (103), the polyhedron is bounded, i.e., $0\leq r_{ij}\leq M$ for some constant $M$. Assume that the polyhedron is nonempty, otherwise the solution is $\infty$. Moreover, a nonempty bounded polyhedron cannot contain a line, thus it must have a basic feasible solution and the optimal solutions are restricted to the corner points. 2. 2.

$\mu_{2}\neq 0$. As $r_{ij}\neq 0$ and $\mu_{1ij}=0$, we have

[TABLE]

Summing (140), (142), (144) and (146) and using (151), we have

[TABLE]

From (154) and $\mu_{2}\neq 0$, we conclude that $c=1$. Moreover, we have $\lambda=\frac{\max(1,\delta)+\mu+\nu\delta}{1+\mu+\nu+\eta}$, where $\mu$, $\nu$ and $\eta$ are defined as:

[TABLE]

Assume that $\eta>0$, therefore $r_{ij}<p_{ij}$ as $p_{ij}\leq\min(q_{i},q_{j})$. On the other hand, $\sum_{i,j}{r}_{ij}=\sum p_{ij}=1$ which contradicts our assumption that $\eta>0$. Thus, $\eta\leq 0$. Define $(\mu^{*},\nu^{*},\eta^{*})$ as follows

[TABLE]

Therefore, from (140), (151) and (154) we conclude Lemma 5 as follows

[TABLE]

D.5 Proof of Lemma 6

In order to prove Lemma 6, consider a leaf node $v_{0}$ and its parent $v_{1}$. From (57), for the node $v_{1}$ we have

[TABLE]

$(\eqref{v1_1}-\eqref{v1_3})$ follow from (57) and the fact that the node $v_{1}$ is niether pruned nor accepted as it is not a leaf node. Therefore, from $\eqref{v1_2}\times\eqref{v1_3}/\eqref{v1_1}$, we conclude that,

[TABLE]

Therefore, from definition of $\Phi(v)$, for the leaf node $v_{0}$ we have

[TABLE]

(167) is true for all the leaf nodes as $(\eqref{v1_1}-\eqref{v1_5})$ were derived independent of the choice of the leaf node.

D.6 Proof of Lemma 7

The proof is straightforward based on proof of Lemma 4. $\sum_{v\in V_{l}}\Phi(v)=1$ is proved by the induction on depth of the tree. For the tree $G$ with $depth(G)=1$, $\sum_{v\in V_{l}}\Phi(v)=1$ is trivial as for the children $v_{ij}$ of the root we have $\Phi(v_{ij})=p_{ij}$ and $\sum_{ij}p_{ij}=1$. Assume that $\sum_{v\in V_{l}}\Phi(v)=1$ is true for all the trees $G$ with $depth(G)\leq depth$. Our goal is to prove that $\sum_{v\in V_{l}}\Phi(v)=1$ is true for all the trees $G$ with $depth(G)=depth+1$. Consider a tree $G$ with $depth(G)=depth+1$ and the tree $G^{\prime}$ obtained by removing all the nodes at depth $depth+1$.

[TABLE]

(168) is a result of definition of $G^{\prime}$, i.e., $G^{\prime}$ is obtained from $G$ by removing all the nodes at depth $depth+1$. (169) is true as for any node $w$ and its children $v_{ij}$ we have $\Phi(w)=\sum_{ij}\Phi(v_{ij})$ which is a result of the fact that $\Phi(v_{ij})=\Phi(w)p_{ij}$. (170) is concluded from the induction assumption, i.e., $\sum_{v\in V_{l}}\Phi(v)=1$ is true for all the trees $G$ with $depth(G)\leq depth$.

D.7 Proof of Lemma 8

For any decision tree which at each node either it is pruned, accepted or branched into $kl$ children, number of nodes in the tree is at most two times number of leaf nodes, i.e., $|V(G)|\leq 2|V_{l}(G)|$. This is true by induction on the depth of the tree. For a tree $G$ with $depth(G)=1$, we have $|V_{l}(G)|=kl$ and $|V(G)|=kl+1$. Therefore, Lemma 8 is true in this case. Assume that $|V(G)|\leq 2|V_{l}(G)|$ is true for all the trees $G$ with $depth(G)\leq depth$. Our goal is to prove that $|V(G)|\leq 2|V_{l}(G)|$ is true for all the trees $G$ with $depth(G)=depth+1$. Consider a tree $G$ with $depth(G)=depth+1$. Consider the tree $G^{\prime}$ obtained by removing all the nodes $v$ where $depth(v)=depth+1$. Assume there are $klr$ of them (each intermediate node has $kl$ children). Therefore, we have $|V(G)|=|V(G^{\prime})|+klr$, $|V_{l}(G)|=|V_{l}(G^{\prime})|+(kl-1)r$ and $|V(G^{\prime})|\leq 2|V_{l}(G^{\prime})|$. This results in $|V(G)|\leq 2|V_{l}(G)|$.

Appendix E Proof of Theorem 4

First of all, note that from $\frac{p_{ij}}{1+\epsilon}\leq p^{\prime}_{ij}\leq{p_{ij}}{(1+\epsilon)}$, we conclude that $\frac{p_{i}^{\mathcal{A}}}{1+\epsilon}\leq{p^{\prime}}_{i}^{\mathcal{A}}\leq{{p}_{i}^{\mathcal{A}}}{(1+\epsilon)}$, $\frac{p_{j}^{\mathcal{B}}}{1+\epsilon}\leq{p^{\prime}}_{j}^{\mathcal{B}}\leq{p_{j}^{\mathcal{B}}}{(1+\epsilon)}$ and $\frac{q_{ij}}{{{(1+\epsilon)}^{2}}}\leq q^{\prime}_{ij}\leq{q_{ij}}{{{(1+\epsilon)}^{2}}}$. Assume that $depth(G)=d$. Define the variables $\Phi(v)$, ${\Psi}^{\mathcal{A}}(v)$, ${\Psi}^{\mathcal{B}}(v)$, $\Psi(v)$, $\alpha(G)$, $\gamma^{\mathcal{A}}(G)$, $\gamma^{\mathcal{B}}(G)$, $\beta(G)$ and $TP$ for the tree with the distribution $p^{\prime}$ as $\Phi^{\prime}(v)$, $\Psi^{\prime\mathcal{A}}(v)$, $\Psi^{\prime\mathcal{B}}(v)$, $\Psi^{\prime}(v)$, $\alpha^{\prime}(G)$, $\gamma^{\prime\mathcal{A}}(G)$, $\gamma^{\prime\mathcal{B}}(G)$, $\beta^{\prime}(G)$ and $TP^{\prime}$. Therefore, from $\eqref{alphaav}-\eqref{gammydef}$ we have

[TABLE]

(172) follows from $\Phi(v)\leq\Phi^{\prime}(v)(1+\epsilon)^{d}$ which is a result of the definition of $\Phi(v)$ in (22), i.e., $\Phi(f(v,a_{i},b_{j}))=\Phi(v)p_{ij},\forall v\in V$. Similarly, we have

[TABLE]

On the other hand, from $\eqref{chagh}-\eqref{approx1}$, we have

[TABLE]

similar to (50), and the inequality $(1-x)^{\frac{c}{x}}<e^{-c}$, the minimum possible value of $\#bands$ to ensure true positive rate $TP$ can be computed as

[TABLE]

Thus, the total complexity is computed as

[TABLE]

where (181) follows from $\eqref{ty_yt1}-\eqref{ty_yt2}$ and the fact that the total complexity for distribution $p^{\prime}$ is $O(N^{\lambda^{*}(p^{\prime})})$. Finally, (182) is obtained from the following lemma in which we prove that the depth of the tree is bounded by $c_{d}\log N$ where $c_{d}$ is a constant depending on the distribution.

Lemma 9

For the decision tree $G$, $depth(G)=c_{d}\log N$ where c_{d}=\min\big{(}\frac{{(\lambda^{*}-1)}}{\log(\max_{i,j}\frac{p_{ij}}{p_{i}^{\mathcal{A}}})},\frac{{(\lambda^{*}-\delta)}}{\log(\max_{i,j}\frac{p_{ij}}{p_{j}^{\mathcal{B}}})}\big{)}.

For proof of Lemma 9, see Appendix E.1.

E.1 Proof of Lemma 9

From (57), for the decision tree $G$ we have

[TABLE]

Our goal here is to prove that for any pruned node $v$ and any accepted node $v$, $depth(v)\leq c_{d}\log N$ where $c_{d}$ is a constant depending on the distribution. Consider a pruned or accepted node $v$. For its parent $w$, we have

[TABLE]

Therefore, we conclude that

[TABLE]

Similarly, we have

[TABLE]

Thus, $c_{d}$ is defined as \min\big{(}\frac{{(\lambda^{*}-1)}}{\log(\max_{i,j}\frac{p_{ij}}{p_{i}^{\mathcal{A}}})},\frac{{(\lambda^{*}-\delta)}}{\log(\max_{i,j}\frac{p_{ij}}{p_{j}^{\mathcal{B}}})}\big{)}. (185) and (187) are true as $p_{0},q_{0}\leq 1$.

Appendix F Pseudo Code

Here, we present the pseudo code for the algorithm for [41] for the case of hamming distance (see Experiment 1).

Appendix G Further Discussion on MIPS

In order to use MIPS to solve this problem, i.e., (2), we need to derive optimal weights $\omega_{ij}$ to minimize the norm $M^{2}$ in [35]. The term $M$ stands for the radius of the space which is computed as follows: M^{2}=\mathbb{E}{\big{(}||x||\big{)}}^{2}+\mathbb{E}{\big{(}||y||\big{)}}^{2}. Therefore, from (69)-(72) we conclude that $M^{2}=\sum_{ij}\left(p_{j}^{\mathcal{B}}\omega_{ij}^{2}+\frac{p_{i}^{\mathcal{A}}{\log^{2}(\frac{p_{ij}}{q_{ij}})}}{\omega_{ij}^{2}}\right)$ which results in optimal \omega_{ij}={\big{(}\frac{p_{i}^{\mathcal{A}}}{p_{j}^{\mathcal{B}}}\big{)}}^{0.25}{\big{(}\mid\log\frac{p_{ij}}{q_{ij}}\mid\big{)}}^{0.5}. On the other hand, for $(x,y)\sim Q(x,y)$ we have

[TABLE]

In order to have nearly one true positive rate and sub-quadratic complexity we need $S_{0}\leq Sd_{KL}(p_{ij}||q_{ij})$ and $cS_{0}\geq-Sd_{KL}(q_{ij}||p_{ij})$ where $d_{KL}$ stands for kullback leibler divergence. Moreover, we should have $M^{2}\geq S\sum_{ij}\sqrt{q_{ij}}|\log(\frac{p_{ij}}{q_{ij}})|$. Setting $c=0$, $S_{0}$ and $M$ as above, the complexity will be more than $1.9$ for any $2\times 2$ probability distribution matrix. The reason is that the transferred data points are nearly orthogonal to each other and this makes it very slow to find maximum inner product using the existing method [35].

Appendix H Complexities of MinHash, LSH-hamming and ForestDSH

In this section, we derive the complexities of MinHash, LSH-hamming and ForestDSH for ${P}_{1}=\begin{bmatrix}0.345&0\\ 0.31&0.345\end{bmatrix}$ for instance. Complexities are computed for the other probability distributions similarly.

H.1 Complexity of MinHash

For MinHash the query complexity is

[TABLE]

where $mh_{1}={\frac{\log\frac{p_{00}}{1-p_{11}}}{\log\frac{q_{00}}{1-q_{11}}}}$, $mh_{2}={\frac{\log\frac{p_{01}}{1-p_{10}}}{\log\frac{q_{01}}{1-q_{10}}}}$, $mh_{3}={\frac{\log\frac{p_{10}}{1-p_{01}}}{\log\frac{q_{10}}{1-q_{01}}}}$ and $mh_{4}={\frac{\log\frac{p_{11}}{1-p_{00}}}{\log\frac{q_{11}}{1-q_{00}}}}$. Similarly for ${P}_{1}$, the per query complexity is derived and is equal to $0.5207$.

H.2 Complexity of LSH-hamming

In the case of LSH-hamming, the query complexity is

[TABLE]

and the storage required for the algorithm is $O(N^{1+\min(\frac{\log(p_{00}+p_{11})}{\log(q_{00}+q_{11})},\frac{\log(p_{01}+p_{10})}{\log(q_{01}+q_{10})})})$. Similarly for ${P}_{1}$, the per query complexity is derived and is equal to $0.4672$.

H.3 Complexity of ForestDSH

From Definition 4, we derive $\lambda^{*}$ as follows

[TABLE]

Note that $\delta=1$ and the per query complexity is equal to $0.4384$.

Appendix I Joint Probability Distributions Learned on Mass Spectrometry Data

The mass spectrometry data for experiment 4, is shown in Figures 10 $(a)$, $(b)$ and $(c)$ in case of $\log Rank$ at base $4$ (a $4\times 4$ matrix), $\log Rank$ at base $2$ (an $8\times 8$ matrix), and no $\log Rank$ transformation (a $51\times 51$ matrix). For the mass spectrometry data shown in Figure 10 $(a)$, the probability distribution $P_{4\times 4}$ is given in (195). Note that, in the case of LSH-hamming the query complexity for these $4\times 4$, $8\times 8$ and $51\times 51$ matrices are $0.901$, $0.890$ and $0.905$, respectively. Similarly, per query complexity for MinHash for these $4\times 4$, $8\times 8$ and $51\times 51$ matrices are $0.4425$, $0.376$ and $0.386$, respectively.

For the mass spectrometry data shown in Figure 10 $a$ and $(b)$, the probability distribution $p(x,y)$ is represented as

[TABLE]

From (56), for $P_{4\times 4}$, $(\mu^{*},\nu^{*},\eta^{*},\lambda^{*})$ are derived as

[TABLE]

Similarly, for $P_{8\times 8}$, we have

[TABLE]

For the mass spectrometry data shown in Figure 10 $(c)$, $(\mu^{*},\nu^{*},\eta^{*},\lambda^{*})$ are derived as

[TABLE]