Sublinear Maximum Inner Product Search using
Concomitants of Extreme Order Statistics

Random projections (RP) refer to the technique of projecting data points in $d$-dimensional spaces onto random $l$-dimensional spaces ($l<d$) via a random matrix $\mathbf{R}\in\mathbb{R}^{l\times d}$. In the reduced $l$-dimensional space, the key data properties, e.g. pairwise Euclidean distances and inner product values, are preserved with high probability. Therefore, we can achieve a high quality approximation answer for MIPS in the reduced dimensional space.

Presenting the point set $\mathbb{X}$ as a matrix $\mathbf{X}\in\mathbb{R}^{d\times n}$, we generate a Gaussian random matrix $\mathbf{R}\in\mathbb{R}^{l\times d}$ whose elements are randomly sampled from the standard normal distribution $N(0,1)$. The signature (i.e. projected representation) of $\mathbf{x}\in\mathbf{X}$ is computed by $\mathbf{R}\mathbf{x}/\sqrt{l}$. Since we study MIPS, we only state the relative distortion bound of inner product values in the reduced space as follows:

The proof of Lemma 1.1 can be found in (Kaban15, ). The constant of this bound on the inner product matches the Johnson-Lindenstrauss lemma bounds on the Euclidean distance (Dasgupta03, ; JL, ), which has been proven to be optimal for linear dimensionality reductions (Larsen17, ).

Locality-sensitive hashing (LSH) (Andoni08, ; Har12, ) is a key algorithmic primitive for similarity search in high dimensions due to the sublinear query time guarantee. Several approaches exploit LSH to obtain sublinear solutions for approximate MIPS (Huang18, ; RangeLSH, ; Shrivastava14, ; SimpleLSH, ). We will review solutions based on SimHash (SimHash, ) since we can use SimHash for both similarity estimation and sublinear search for MIPS.

SimHash. Given a Gaussian random vector $\mathbf{r}\in\mathbb{R}^{d}$ whose elements are randomly drawn from the $N(0,1)$, a SimHash function of $\mathbf{x}$ is $h(\mathbf{x})=\operatorname{sgn}\left(\mathbf{x}^{\top}\mathbf{r}\right)$. Denote by $0\leq\theta_{\mathbf{x}\mathbf{q}}\leq\pi$ the angle between two vectors $\mathbf{x}$ and $\mathbf{q}$, the seminal work of Goemans and Williamson (Goemans, ) and Charikar (SimHash, ) show that

Since $\arccos()$ is a monotonically decreasing function of $\mathbf{x}^{\top}\mathbf{q}$ when $\|\mathbf{x}\|=\|\mathbf{q}\|=1$, SimHash is a LSH family for the inner product similarity on a unit sphere. Exploiting the binary representation of SimHash values, we can efficiently estimate the inner product $\mathbf{x}^{\top}\mathbf{q}$ with the fast Hamming distance computation using built-in functions of compilers. We note that SimHash binary code can be seen as a quantization of Gaussian RP since it keeps the sign of each projected value.

SimHash-based solutions for MIPS. Since the change of $\|\mathbf{q}\|$ does not affect the result of MIPS, we can assume $\|\mathbf{q}\|=1$ without loss of generality. Lemma 1.2 shows that the hash collision depends on both $\mathbf{x}^{\top}\mathbf{q}$ and $\|\mathbf{x}\|$. Therefore, by storing the 2-norm $\|\mathbf{x}\|$, we can use SimHash for estimating $\mathbf{x}^{\top}\mathbf{q}$. However, the dependency on such 2-norms makes MIPS more challenging to guarantee a sublinear query time. Since inner product is not a metric, SimHash-based solutions have to convert MIPS to the nearest neighbor search by applying order preserving transformations to ensure that both data and query are on a unit sphere.

SimpleLSH (SimpleLSH, ) proposes asymmetric transformations : $\mathbf{q}\mapsto\mathbf{q}^{\prime}=\{\mathbf{q},0\}$ and $\mathbf{x}\mapsto\mathbf{x}^{\prime}=\{\mathbf{x}/M,\sqrt{1-\|\mathbf{x}\|^{2}/M^{2}}\}$ where $M$ is the maximum 2-norm of all points in $\mathbf{X}$. Since the inner product order is preserved and $\|\mathbf{x}^{\prime}\|=\|\mathbf{q}^{\prime}\|=1$, we can exploit SimHash for both similarity estimation and sublinear search for MIPS.

It is clear that the inner product values in the transformation space will be scaled down by a factor of $M$. Furthermore, the top-$k$ inner product values are often very small compared to the vector norms in high dimensions. This means that the MIPS values and the distance gaps between “close” and “far apart” points in the transformation space are arbitrary small. Therefore, we need to use significantly large code length to achieve a reasonable search recall. In addition, the standard $(l,L)$-parameterized (where $l$ is the number of concatenating hash functions and $L$ is number of hash tables) bucketing algorithm (Andoni08, ) demands a huge space usage (i.e. large $L$) to guarantee a sublinear query time. In other words, the LSH performance will be degraded in the transformation space.

RangeLSH (RangeLSH, ) handles this problem by partitioning the data into several partitions and applying SimpleLSH on each partition. While such distance gaps between “close” and “far apart” points in the LSH framework are slightly improved, the number of hash tables $L$ needs to be scaled proportional to the number of partitions to achieve a reasonable search recall. Nevertheless, the subquadratic time and space cost of building LSH tables in order to guarantee a sublinear query time will be a bottleneck on many big data applications.

The paper introduces a specific dimensionality reduction, called CEOs, based the theory of concomitants of extreme order statistics. While we also construct signatures for both data and query points using $D$ Gaussian random projections, we only use a small subset of $D$ dimensions to estimate inner products. This subset of size $s\ll D$ corresponds to random projections associated with the $s$ extreme values (i.e. maximum or minimum) of the query signature.

The geometric intuition is that among $D$ Gaussian random projection vectors $\mathbf{r}_{i},1\leq i\leq D$, the random projection associated with $\mathbf{r}_{(1)}$, the closest one to the query $\mathbf{q}$, will preserve the best inner product order. Since $\mathbf{r}_{(1)}$ is the closest vector to $\mathbf{q}$, the value $\mathbf{q}^{\top}\mathbf{r}_{(1)}$ is maximum among $D$ signature values $\mathbf{q}^{\top}\mathbf{r}_{i}$. Figure 1 shows a high level illustration of how CEOs works. We note that ones can also use the random projection associated with $\mathbf{r}_{(D)}$, the furthest one to $\mathbf{q}$ due to the symmetry of Gaussian distribution.

The technical difficulty is to provide a good estimator for $\mathbf{x}^{\top}\mathbf{q}$ given a small subset of projected values associated with the $s$ random projection vectors closest and furthest to $\mathbf{q}$. We observe that this subset of size $s\ll D$ is the concomitants of extreme $s$th normal order statistics associated with the query signature. Leveraging the theory of concomitants of extreme order statistics, we provide a surprisingly simple and asymptotically unbiased estimate for $\mathbf{x}^{\top}\mathbf{q}$. That yields very fast and accurate MIPS solvers.

The key feature of CEOs is that we only need to use significantly small $s$ projected dimensions (often $s=10$ in our benchmark) among $D$ signature dimensions. Most importantly, in contrast to traditional dimensionality reductions, CEOs only uses the sign of a small subset of query signature (i.e. projections associated with maximum and minimum values) to estimate inner products. This means that we can precompute and rank inner product estimators of all data before querying to significantly reduce the query time. We show that this cost of constructing indexes matches the lower bound of the $(1+\epsilon)$-approximate MIPS in unit sphere (Andoni06, ).

Under a mild condition, our theoretical analysis shows that we can achieve a sublinear query time for MIPS with search recall guarantees given an ${\mathcal{O}}\left((D/s)^{s}\right)$ indexing space complexity. To deal with such exponential time and space complexity, we introduce several CEOs variants which require $\tilde{{\mathcal{O}}}\left(dn\right)$ time and space for indexing and can answer top-$k$ MIPS with very high search recalls. A summary of our contributions is as follows:

1. (1)

   We propose a novel dimensionality reduction method based on the theory of concomitants of extreme order statistics, named CEOs for top-$k$ MIPS. We theoretically and empirically show that CEOs provides better estimate accuracy and hence higher top-$k$ MIPS recall than SimHash (SimHash, ) and Gaussian random projections (Kaban15, ).

2. (2)

   Inheriting from the asymptotic properties of concomitants of extreme $s$th order statistics, CEOs yields a sublinear query time given ${\mathcal{O}}\left((D/s)^{s}\right)$ time and space complexity for building the index. The search recall can be theoretically guaranteed under a mild condition, and the exponential space usage is optimal for MIPS on a unit sphere.

3. (3)

   We propose practical CEOs variants, including sCEOs-TA and coCEOs, that build the index in $\tilde{{\mathcal{O}}}\left(dn\right)$ time. sCEOs-TA exploits the threshold algorithm (TA, ) for solving MIPS. coCEOs is a data and dimension co-reduction technique and often outperforms sCEOs-TA on high search recall requirements. Especially, coCEOs can be viewed as a budgeted MIPS solver with a parameter $B=o(n)$ that answers top-$k$ MIPS in $o(n)$ time.

4. (4)

   Our proposed algorithms are very simple to implement and require few lines of codes. Empirically, on the inferior choice that uses the concomitants of the maximum order statistics, CEOs requires sublinear index space and outperforms LSH bucket algorithms regarding both efficiency and accuracy. Both sCEOs-TA and coCEOs outperform competitive MIPS solvers on answering top-10 MIPS and achieve at least 100x speedup compared to the bruteforce search with the accuracy at least 90% on many real-world large-scale data sets.

Let $(U,V)$ be bivariate normal $N\left(\mu_{U},\mu_{V},\sigma_{U}^{2},\sigma_{V}^{2},\rho\right)$ where $\rho$ is the correlation coefficient, and $\mu_{U},\mu_{V},\sigma_{U}^{2},\sigma_{V}^{2}$ are the means and variances of $U$ and $V$, respectively. Let $(U_{i},V_{i})$ be a bivariate normal sample, we can write

where $U_{i}$ and $\epsilon_{i}$ are mutually independent. Furthermore, $\epsilon_{i}$ is normal with $\textrm{\bf E}\left[\epsilon_{i}\right]=0,\textrm{\bf Var}\left[\epsilon_{i}\right]=\sigma_{V}^{2}(1-\rho^{2})$.

(Bivariate_Book, , Section 4.7) shows that the conditional distribution of $V$ given $U=x$ is normal with the mean and variance as follows.

Let $(U_{1},V_{1}),(U_{2},V_{2}),\ldots,(U_{D},V_{D})$ be $D$ random samples from a bivariate normal distribution $N\left(\mu_{U},\mu_{V},\sigma_{U}^{2},\sigma_{V}^{2},\rho\right)$. We order these samples based on the $U$-value. Given the $r$th order statistic $U_{(r)}$, the $V$-value associated with $U_{(r)}$ is called concomitant of the $r$th order statistic and denoted by $V_{[r]}$. Let $\epsilon_{[r]}$ denote the specific $\epsilon_{i}$ associated with $U_{(r)}$, we have the following equation:

The seminal work of David and Galambos (CEO_Paper, ) establishes the following properties of concomitants of normal order statistics.

In this work, we are interested in the concomitants of a few $r$th and $(D-r)$th order statistics where $r$ is small, for example $V_{[1]}$ and $V_{[D]}$. The asymptotic distribution of these random variables with a sufficiently large $D$ has been studied in statistics during the last decades (CEO_Book, ; CEO_Paper, ) with a considerable use in survival analysis.

Let $X_{[1]}$ be the concomitant of the first (maximum) order statistic $Q_{(1)}$. Applying Lemma 2.2, we have the following properties of $X_{[1]}$:

We note that $Q_{(1)}$ is the largest variable among $D$ independent standard normal variables. When $D$ is sufficiently large, (maxGauss_Book, , Chapter 8) and (maxGauss, ) show that $\textrm{\bf E}\left[Q_{(1)}\right]\rightarrow\sqrt{2\log{D}}$ and $\textrm{\bf Var}\left[Q_{(1)}\right]\rightarrow 0$, and the rate of convergence is ${\mathcal{O}}\left(1/\log{D}\right)$. Using $\xrightarrow{D}$ as a proxy for asymptotic results with a sufficiently large $D$, we have

In order to use the concomitant of the first order statistics for estimating inner product value, we define $Z_{ceo}=X_{[1]}/\sqrt{2\log{D}}$. It is straightforward that $\textrm{\bf E}\left[Z_{ceo}\right]\xrightarrow{D}\mathbf{x}^{\top}\mathbf{q}$ and

We note that the behavior of $X_{[1]}$ is consistent with the result of Lemma 2.1 where $x$ are the extreme values of standard normal distribution, e.g. around $\sqrt{2\log D}$. While the variance is stable, the expectation is scaled up by a factor of $\sqrt{2\log D}$, which yields highly accurate estimates for MIPS. For notational simplicity, we name the method using $X_{[1]}$ for estimating inner products as CEOs.

Comparison with LSH: We now compare the inner product estimators provided by SimHash and CEOs. Define $Z=1$ if $h(\mathbf{x})=h(\mathbf{q})$; otherwise 0. By the LSH property, we have $\textrm{\bf E}\left[Z\right]=1-\theta_{\mathbf{x}\mathbf{q}}/\pi$ and $\textrm{\bf Var}\left[Z\right]=\left(\theta_{\mathbf{x}\mathbf{q}}/\pi\right)\left(1-\theta_{\mathbf{x}\mathbf{q}}/\pi\right)$. Using the Taylor series with $\arccos{x}\approx\pi/2-x$ for $-1\leq x\leq 1$, we have

Define $Z_{lsh}=(Z-1/2)\pi\|\mathbf{x}\|$, it is clear that $\textrm{\bf E}\left[Z_{lsh}\right]\approx\mathbf{x}^{\top}\mathbf{q}$ and

Equations 1 and 2 show that SimHash needs to use approximately $\frac{\pi^{2}}{2}\log{D}\approx 5\log{D}$ bits to achieve a similar accuracy as CEOs with a sufficiently large $D$. As elaborated in Section 1.2, SimpleLSH applies SimHash for MIPS by scaling down all inner products with a normalization factor $M$ where $M$ is the maximum 2-norm of the data. This means that SimpleLSH needs to use approximately $5M^{2}\log{D}$ bits to have a similar accuracy as CEOs.

Comparison with Gaussian RP: We now bound the estimate error provided by $X_{[1]}$ and compare with the concentration bound provided by the Gaussian RP. We investigate the asymptotic behavior of concomitants of extreme order statistics. The theory of extreme order statistics (CEO_Paper, ) states that the distribution of the extreme case $X_{[1]}$ converges in probability to the normal distribution $N\left(\mathbf{x}^{\top}\mathbf{q}\sqrt{2\log{D}}\,,\|\mathbf{x}\|^{2}-\left(\mathbf{x}^{\top}\mathbf{q}\right)^{2}\right)$ when $D\rightarrow\infty$. Recall that $Z_{ceo}=X_{[1]}/\sqrt{2\log{D}}$ and hence $Z_{ceo}\xrightarrow{D}N\left(\mathbf{x}^{\top}\mathbf{q}\,,\frac{\|\mathbf{x}\|^{2}-\left(\mathbf{x}^{\top}\mathbf{q}\right)^{2}}{2\log{D}}\right)\,$.

We will use the following Chernoff bounds (Chernoff, ) to bound the relative estimate error induced by $X_{[1]}$.

Applying the Chernoff bounds for the random variable $Z_{ceo}$, given any $\epsilon\geq 0$ we have:

When $D$ is sufficiently large, the error probability will be arbitrarily small, and hence the estimate will be very accurate. Using the asymptotic distribution of $X_{[1]}$, our concentration bound has the similar form of Lemma 1.1 with $l=8\log{D}/\sin^{2}{(\theta_{\mathbf{x}\mathbf{q}})}$ projections. While in the worst case where $\theta_{\mathbf{x}\mathbf{q}}\rightarrow\pi/2$, Gaussian RP shares the similar bound as CEOs with ${\mathcal{O}}\left(\log{D}\right)$. However for the case of $\theta_{\mathbf{x}\mathbf{q}}\rightarrow 0$, CEOs offers a much tighter bound than Gaussian RP due to $\sin^{2}{(\theta_{\mathbf{x}\mathbf{q}})}\rightarrow 0$.

MIPS algorithms: In general, compared to CEOs, SimHash and Gaussian RP can achieve similar performance for MIPS by using ${\mathcal{O}}\left(\log{D}\right)$ bits and ${\mathcal{O}}\left(\log{D}\right)$ random projections for constructing signatures, respectively. However, after constructing signatures, these approaches have to perform ${\mathcal{O}}\left(n\right)$ inner product estimation, which is clearly a computational bottleneck for MIPS. In contrast, CEOs does not need this costly estimation step since we can sort all the points based on their projection values in advance. At the query phase, after executing RP to compute the dimension index of $Q_{(1)}$, we can simply return $k$ point indexes corresponding to the top-$k$ concomitants associated with $Q_{(1)}$ as an approximate top-$k$ MIPS. This key property makes CEOs more efficient than both Gaussian RP and LSH-based solutions on answering approximate MIPS.

Theoretically, the theory of concomitants of extreme order statistics requires $D\rightarrow\infty$. However, we observe in practice that $D=2^{10}$ suffices for many real-world data sets. We present a brief experiment to compare the MIPS accuracy using $X_{[1]}$ called CEOs with $D=2^{10}$, SimHash, SimpleLSH and Gaussian RP. We use the precision measure $P@k$ which computes the search recall of the top-$k$ answers provided by the estimation algorithms. Figure 2 shows the average measure $P@10$ and $P@20$ for $k=10$ and $k=20$ of these 4 methods on Nuswide over 100 queries. The result is consistent with our analysis. CEOs with $D=2^{10}$ outperforms Gaussian RP with $l=10$, and achieves higher accuracy than SimpleLSH with $l\leq 128$ and SimHash with $l<32$ bits code.

We observe that the Gaussian distribution is symmetric. This means that we can use both $X_{[1]}$ and $X_{[D]}$ corresponding to the maximum $Q_{(1)}$ and minimum $Q_{(D)}$ order statistics for estimating $\mathbf{x}^{\top}\mathbf{q}$. One natural question is: “Can we use more concomitants and are they independent so that we can boost the accuracy with the standard Chernoff bounds?”.

Indeed, we have a positive answer by investigating the asymptotic independence of concomitants of extreme order statistics. Let $s_{0}$ be a fixed positive integer and $D\rightarrow\infty$, the seminal work of David and Galambos (CEO_Paper, ) shows that the concomitants $X_{[r]}$ where $1\leq r\leq s_{0}$ are asymptotically independent and normal as follows:

By the symmetry of Gaussian distribution, we can use the concomitants $X_{[r]}$ and $X_{[D-r+1]}$ where $1\leq r\leq s_{0}$ to boost the accuracy of the estimator of $\mathbf{x}^{\top}\mathbf{q}$. In other words, ones can view these $s=2s_{0}$ concomitants as a specific dimension reduction since the sum of these values provides unbiased estimates of inner products.

By using the average of these $s$ concomitants to estimate $\mathbf{x}^{\top}\mathbf{q}$, we achieve an asymptotic unbiased estimator with the variance decreased by a factor of $1/s$. Hence, the concentration bound is tighter by a factor of $1/\exp(s)$.

Figure 3 (a) confirms our results on Nuswide by showing the average precision $P@10$ of CEOs on top-10 MIPS when varying $s$ and $D$. When fixing $s$ and $D$ is sufficiently large, e.g. $\geq 512$, the precision marginally increases. However, increasing $s$ yields a substantial rise of precision. Especially, when $D=512$ and $s=4$, $P@10$ is nearly 90%, which is higher than that of SimHash with 128 bits code shown in Figure 2 (a).

This subsections discusses the optimality of CEOs on answering the $(1+\epsilon)$-approximate MIPS on a unit sphere given a small $\epsilon$. Given a point set $\mathbb{X}\subset\mathbb{R}^{d}$ of size $n$ and a query point $\mathbf{q}\in\mathbb{R}^{d}$ on a unit sphere, it is clear that $\operatorname*{arg\,max}_{\mathbf{x}\in\mathbb{X}}{\mathbf{x}^{\top}\mathbf{q}}=\operatorname*{arg\,min}_{\mathbf{x}\in\mathbb{X}}{\|\mathbf{x}-\mathbf{q}\|}$. Therefore, we investigate the result of theorem 4.1 for a decision version of the approximate nearest neighbor search over the Euclidean space on a unit sphere. The decision version is to build a data structure that answers $(1+\epsilon)$-approximate near neighbor search as follows.

• •

  If there is $\mathbf{x}\in\mathbb{X}$ such that $\|\mathbf{x}-\mathbf{q}\|\leq r$, answer YES.

• •

  If there is no $\mathbf{x}\in\mathbb{X}$ such that $\|\mathbf{x}-\mathbf{q}\|\leq(1+\epsilon)r$, answer NO.

Andoni et al. (Andoni06, ) analyze the lower bound of the decision problem given a constant query time (measured by the number of probes to the data structure) on the Euclidean space. In particular, any algorithm uses a constant number of probes to the data structure must use $n^{\Omega(1/\epsilon^{2})}$ space. Since a high search recall demands very small $\epsilon$, this result presents a fundamental difficulty for any practical data structures that can answer the nearest neighbor search very accurate and fast.

An upper bound for $(1+\epsilon)$-nearest neighbor search over the Hamming space is proposed in (Kushilevitz00, ) that use $n^{{\mathcal{O}}(1/\epsilon^{2})}$. This work complements that result by showing a data structure which achieves a constant query time and uses $n^{{\mathcal{O}}(1/\epsilon^{2})}$ over the Euclidean space on a unit sphere. We show it by deriving and bounding the value $\alpha_{\mathbf{y}}$ on Theorem 4.1 where $\|\mathbf{x}-\mathbf{q}\|\leq r$ and $\|\mathbf{y}-\mathbf{q}\|\geq(1+\epsilon)r$ as follows.

Hence,

On a unit sphere and for a fixed distance $0\leq r\leq 2$, we have $\alpha_{\mathbf{y}}=\Theta(\epsilon)$. From Equation 4, we have to use $D=n^{1/\alpha_{\mathbf{y}}^{2}}=n^{\Theta(1/\epsilon^{2})}$ to answer the decision version of $(1+\epsilon)$-approximate near neighbor search with high probability. This means that CEOs can precompute the results of $(1+\epsilon)$-approximate near neighbor search for all queries using $n^{\Theta(1/\epsilon^{2})}$ space and answer the query with a constant probes. This matches the lower bound on the space usage provided by (Andoni06, ).

Consider that we have a small budget of $b$ inner product computations for post-processing to achieve higher MIPS accuracy. For each projected dimension, we only need to maintain top-$b$ largest concomitants among $n$ concomitants corresponding to $n$ points. For querying, we only use these $b$ concomitants associated with $Q_{(1)}$. Since any MIPS solver can use post-processing, we will consider $b={\mathcal{O}}\left(1\right)$ to simplify the analysis. Algorithm 1 shows how CEOs with concomitants of $Q_{(1)}$, named 1CEOs, works.

Complexity: It is clear that building the index takes ${\mathcal{O}}\left(dDn\right)$ time. The index space is only ${\mathcal{O}}\left(dD\right)$ since each of $D$ dimensions stores $b$ points, and the query time is ${\mathcal{O}}\left(dD\right)$. When $D=o(n)$, 1CEOs answers MIPS in both sublinear space and time.

Error analysis: Equation 3 indicates that, given a sufficiently large $D$, $X_{[1]}$ corresponding to the top-$k$ MIPS tends to be ranked at the top positions on the list $L_{j_{(1)}}$ corresponding to $Q_{(1)}$. We note that by simply increasing $b$ inner products in post-processing, we can achieve higher accuracy of top-$k$ MIPS since we allow larger gap $Y_{[1]}-X_{[1]}$. While we cannot theoretically guarantee the performance of 1CEOs, our empirical results show that 1CEOs outperforms the sublinear LSH bucket algorithm regarding both space and time of indexing and querying for MIPS.

Since we can use $s=2s_{0}$ concomitants associated with $Q_{(r)}$ and $Q_{(D+1-r)}$ where $1\leq r\leq s_{0}$ to boost the estimate accuracy, Theorem 4.1 shows that we can have a sublinear MIPS with guarantees.

Since the estimation does not use the values of query signatures (except the dimension order), we can precompute and sort all possible estimates before querying. For example, in order to use $X_{[1]}$ and $X_{[D]}$ for estimating $\mathbf{x}^{\top}\mathbf{q}$, after executing $D$ random projections, we can precompute the difference of all pairs of dimensions and sort the data based on these difference values. We compute the difference $X_{[1]}-X_{[D]}$ because of the symmetry of normal bivariate distribution. Algorithm 2, named sCEOs, generalizes this observation to exploit concomitants of extreme $s$th order statistics for MIPS.

Complexity and error analysis: Building the index takes ${\mathcal{O}}\left(\left(D/s\right)^{s}n\right)$ time and ${\mathcal{O}}\left(\left(D/s\right)^{s}+dn\right)$ space. This is because we have $\binom{D}{s_{0}}$ different sets $I$ and each of $I$ has $\binom{D-s_{0}}{s_{0}}$ different sets $J$. sCEOs has ${\mathcal{O}}\left(dD\right)$ query time as similar as 1CEOs. Setting $D=n^{1/c}$ for any $c>1$, Theorem 4.1 states that sCEOs can have a sublinear query time given that the assumption $\alpha^{2}_{\mathbf{y}}>c/s_{0}$ holds for all $\mathbf{y}$. Furthermore, sCEOs returns exact top-$k$ MIPS with probability at least $1-1/n$. Figure 3 (a) demonstrates that with $s=4$, $D=512$ and $b=10$, sCEOs can achieve nearly 90% search recall for top-10 MIPS on Nuswide even without post-processing.

Performance simulation: One of the nice features of sCEOs is that we can compute the sCEOs search recall by implementing it as a specific dimensionality reduction. We name this approach as sCEOs-Est since it estimates $n$ inner products by the sum of concomitants of extreme $s$th order statistics of the query signature. While sCEOs-Est estimates $n$ inner products in ${\mathcal{O}}\left(sn\right)$ time, it still runs significantly faster than both Gaussian RP and bruteforce search since addition operators are often much faster than multiply-add operators. More important, we can use sCEOs-Est to simulate the sublinear search performance of sCEOs before implementing it. This property makes sCEOs very useful in practice while LSH bucket algorithms do not have.

While the querying time is sublinear, the space and time complexity for building the index of sCEOs blow up by a factor of $(D/s)^{s}$, which will be the computational bottleneck of many big data applications. Therefore, we might need to use sCEOs-Est for large data sets.

We observe that sCEOs-Est estimates inner product values by the sum of $s$ positive concomitants and extracts top-$b$ indexes with the largest estimates for post-processing. Therefore, ones can speed up sCEOs-Est by using the well-known threshold algorithm (TA) (TA, ). Algorithm 3, named sCEOs-TA, shows how we can exploit the TA algorithm for speeding up sCEOs-Est.