Index-based, High-dimensional, Cosine Threshold Querying with Optimality Guarantees

TL;DR

This paper introduces optimized algorithms for cosine threshold queries on high-dimensional data, leveraging data normalization and skewness to improve efficiency with strong optimality guarantees.

Contribution

It presents novel traversal strategies and stopping conditions for index-based cosine similarity queries, achieving near-optimal performance under realistic data assumptions.

Findings

Improved stopping condition for index traversal using vector normalization

A novel traversal strategy exploiting data skewness

Strong near-optimal performance guarantees under skewness assumptions

Abstract

Given a database of vectors, a cosine threshold query returns all vectors in the database having cosine similarity to a query vector above a given threshold {\theta}. These queries arise naturally in many applications, such as document retrieval, image search, and mass spectrometry. The present paper considers the efficient evaluation of such queries, providing novel optimality guarantees and exhibiting good performance on real datasets. We take as a starting point Fagin's well-known Threshold Algorithm (TA), which can be used to answer cosine threshold queries as follows: an inverted index is first built from the database vectors during pre-processing; at query time, the algorithm traverses the index partially to gather a set of candidate vectors to be later verified for {\theta}-similarity. However, directly applying TA in its raw form misses significant optimization opportunities. Indeed, we first show that one can take advantage of the fact that the vectors can be assumed to be normalized, to obtain an improved, tight stopping condition for index traversal and to efficiently compute it incrementally. Then we show that one can take advantage of data skewness to obtain better traversal strategies. In particular, we show a novel traversal strategy that exploits a common data skewness condition which holds in multiple domains including mass spectrometry, documents, and image databases. We show that under the skewness assumption, the new traversal strategy has a strong, near-optimal performance guarantee. The techniques developed in the paper are quite general since they can be applied to a large class of similarity functions beyond cosine.

Figures

Tables

Table 1. Table 1: Summary of theoretical results for the near-convex case.

	Stopping Condition		Traversal Strategy
	Baseline	This work	Baseline	This work
Inner Product	Tight		$m\cdot \mathrm{𝖮𝖯𝖳}$	$\mathrm{𝖮𝖯𝖳}+c$
Cosine	Not tight	Tight	$\mathrm{𝖭𝖠}$	$\mathrm{𝖮𝖯𝖳}\left(\theta -ϵ\right)+c$

Equations

$$(a)\sum _{i=1}^d{s}_i\cdot q_i\ge \theta ,(b)s_i\le L_i[b_i],\forall i\in [d],\text{ and}(c)\sum _{i=1}^d{s}_i^2=1.\text{vspace}-1.5mm$$

$$\mathsf{MS}(L[\mathbf{b}])=\sum _{i=1}^d\min \{q_i\cdot \tau ,L_i[b_i]\}\cdot q_i$$

$$\vspace{-1.5mm}\begin{array}[]{ll@{}ll@{}ll@{}ll}\text{maximize}&\displaystyle\sum_{i=1}^{d}s_{i}\cdot q_{i}&\text{subject to}&\displaystyle\sum_{i=1}^{d}s_{i}^{2}=1\quad\text{ and }\quad s_{i}\leq L_{i}[b_{i}],\quad\forall i\in[d]\end{array}\vspace{-1.5mm}$$

$$L(\mathbf{s},\mu ,\lambda )=\sum _{i=1}^d{s}_i{q}_i-\sum _{i=1}^d{\mu }_i(L_i[b_i]-s_i)-\lambda \left(\sum _{i=1}^d{s}_i^2-1\right)$$

$$\begin{array}[]{ll@{}ll}\nabla_{\mathbf{s}}L(\mathbf{s},\mu,\lambda)=0&\text{(Stationarity)}\\ \mu_{i}\geq 0,\ \forall\ i\in[d]&\text{(Dual feasibility)}\\ \mu_{i}(L_{i}[b_{i}]-s_{i})=0,\ \forall\ i\in[d]&\text{(Complementary slackness)}\end{array}$$

$${\mathcal{T}}_{\mathsf{MR}}(\mathbf{b})=\underset{1\le i\le d}{\mathrm{argmax}}\left(F(L[\mathbf{b}])-F(L[\mathbf{b}+\text{mathds}{1}_i])\right)=\underset{1\le i\le d}{\mathrm{argmax}}\left(f_i(L_i[b_i])-f_i(L_i[b_i+1])\right)$$

$${\mathcal{T}}_{\mathsf{HL}}(\mathbf{b})=\underset{1\le i\le d}{\mathrm{argmax}}({\tilde{\Delta }}_i[b_i]).$$

$$F(L[{\mathbf{b}}^{\ast }])\ge \tilde{F}(L[{\mathbf{b}}^{\ast }])-{ϵ}_1\ge \tilde{F}(L[{\mathbf{b}}_l])-{ϵ}_1\ge F(L[{\mathbf{b}}_l])-{ϵ}_1-{ϵ}_2\ge \theta -{ϵ}_1-{ϵ}_2.\text{vspace}-1mm$$

$$\tilde{F}(L[\mathbf{b}])=\sum _{i=1}^d\min \{q_i\cdot \tilde{\tau },L_i[b_i]\}\cdot q_i\text{vspace}-1mm$$

$$ϵ\le \max \{0,\tilde{\tau }-1/\mathsf{MS}(L[{\mathbf{b}}_l])\}+\mathsf{MS}(L[{\mathbf{b}}_l])-\tilde{F}(L[{\mathbf{b}}_l]).\text{vspace}-1mm$$

$$\big{(}q_{i}\cdot\tilde{\tau}-L_{i}[j_{k}]\big{)}\ /\ j_{k}\geq\big{(}L_{i}[j_{k}]-L_{i}[j_{k+1}]\big{)}\ /\ (j_{k+1}-j_{k}).\vspace{-1mm}$$

$${\phi }_{\mathsf{BL}}(\mathbf{b})=\left(\sum _{i=1}^d{q}_i\cdot L_i[b_i]<\theta \right).$$

$${\mathcal{T}}_{\mathsf{BL}}(\mathbf{b})=\mathsf{nz}[|\mathbf{b}|\mathsf{mod}m],$$

$$ub(\tilde{\mathbf{s}})=\sum _{\tilde{\mathbf{s}}[i]\ne \perp }\tilde{\mathbf{s}}[i]\cdot \mathbf{q}[i]+\sqrt{1-\sum _{\tilde{\mathbf{s}}[i]\ne \perp }\tilde{\mathbf{s}}[i{]}^2}\cdot \sqrt{1-\sum _{\tilde{\mathbf{s}}[i]\ne \perp }\mathbf{q}[i{]}^2},$$

$$lb(\tilde{\mathbf{s}})=\sum _{\tilde{\mathbf{s}}[i]\ne \perp }\tilde{\mathbf{s}}[i]\cdot \mathbf{q}[i]+\sqrt{1-\sum _{\tilde{\mathbf{s}}[i]\ne \perp }\tilde{\mathbf{s}}[i{]}^2}\cdot \underset{\tilde{\mathbf{s}}[i]=\perp }{\min }\mathbf{q}[i].$$

$$lb$$

$$ub$$

$$\sum _{i=1}^k\mathbf{s}[i]\cdot \mathbf{q}[i]+\sum _{i=k+1}^d\mathbf{s}[i]\cdot \mathbf{q}[i]\ge \theta +\sqrt{1-c},$$

$$\sum _{i=k+1}^d\mathbf{s}[i]\cdot \mathbf{q}[i]$$

$$=\sqrt{1-\sum _{i=1}^k\mathbf{s}[i{]}^2}\cdot \sqrt{1-\sum _{i=1}^k\mathbf{q}[i{]}^2}$$

$$\le \sqrt{1-c}.$$

$$\sum _{i=1}^k\mathbf{s}[i]\cdot \mathbf{q}[i]+\sqrt{1-\sum _{i=1}^k\mathbf{s}[i{]}^2}\cdot \sqrt{1-\sum _{i=1}^k\mathbf{q}[i{]}^2}\le \theta .$$

$${\phi }_{\mathsf{BL}}(\mathbf{b})=\left(\sum _{i=1}^d{q}_i\cdot L_i[b_i]<\theta \right)$$

$$s_i=\{\begin{array}{ll}{L}_i[b_i],& \tau \ge \frac{L_i[b_i]}{q_i};\\ q_i\cdot \tau ,& \text{otherwise}.\end{array}$$

$$\mathsf{eval}(k,\tau )=\sum _{i=1}^d{s}_i^2=\sum _{i=1}^k{L}_i[b_i{]}^2+\sum _{i=k+1}^d{q}_i^2\cdot {\tau }^2,$$

$$\sum _{i=1}^k{L}_i[b_i{]}^2+\sum _{i=k+1}^d{q}_i^2\cdot {\tau }^2=1\Rightarrow \tau ={\left(\frac{1-{\sum }_{i=1}^k{L}_i[b_i{]}^2}{1-{\sum }_{i=1}^k{q}_i^2}\right)}^{1/2}.$$

$$\mathsf{MS}(L[\mathbf{b}])=\sum _{i=1}^k{L}_i[b_i]\cdot q_i+(1-\sum _{i=1}^k{q}_i^2)\cdot \tau .$$

$$s_i\ge L_i[a_i]\vee s_i\le L_i[b_i],$$

$$(s_i-L_i[a_i])\cdot (s_i-L_i[b_i])\ge 0$$

$$\begin{array}[]{ll@{}ll@{}ll}\text{maximize}&\displaystyle\sum_{i=1}^{d}s_{i}\cdot q_{i}&\\ \text{subject to}&\displaystyle\sum_{i=1}^{d}s_{i}^{2}=1,&\\ &(s_{i}-L_{i}[a_{i}])\cdot(s_{i}-L_{i}[b_{i}])\geq 0,&\quad\text{ for }i\in[d].\end{array}$$

$$F(\mathbf{s},\mu ,\lambda )=\sum _{i=1}^d{s}_i{q}_i+\sum _{i=1}^d{\mu }_i(s_i-L_i[a_i])\cdot (s_i-L_i[b_i])-\lambda \left(\sum _{i=1}^d{s}_i^2-1\right)$$

$$\begin{array}[]{ll@{}ll}\nabla_{\mathbf{s}}F(\mathbf{s},\mu,\lambda)=0&\text{(Stationarity)}\\ \mu_{i}\geq 0,\ \forall\ i\in[d]&\text{(Dual feasibility)}\\ \mu_{i}(s_{i}-L_{i}[a_{i}])(s_{i}-L_{i}[b_{i}])=0,\ \forall\ i\in[d]&\text{(Complementary slackness)}\end{array}$$

Full text

1 Introduction

Given a database of vectors, a cosine threshold query asks for all database vectors with cosine similarity to a query vector above a given threshold.

This problem arises in many applications including document retrieval [14], image search [46], recommender systems [51] and mass spectrometry. For example, in mass spectrometry, billions of spectra are generated for the purpose of protein analysis [1, 47, 71]. Each spectrum is a collection of key-value pairs where the key is the mass-to-charge ratio of an ion contained in the protein and the value is the intensity of the ion. Essentially, each spectrum is a high-dimensional, non-negative and sparse vector with $\sim$2000 dimensions where $\sim$100 coordinates are non-zero.

Cosine threshold queries play an important role in analyzing such spectra repositories. Example questions include “is the given spectrum similar to any spectrum in the database?”, spectrum identification (matching query spectra against reference spectra), or clustering (matching pairs of unidentified spectra) or metadata queries (searching for public datasets containing matching spectra, even if obtained from different types of samples). For such applications with a large vector database, it is critically important to process cosine threshold queries efficiently – this is the fundamental topic addressed in this paper.

Definition 1.1** (Cosine Threshold Query).**

Let $\mathcal{D}$ be a collection of high-dimensional, non-negative vectors; $\mathbf{q}$ be a query vector; $\theta$ be a threshold $0<\theta\leq 1$. Then the cosine threshold query returns the vector set $\mathcal{R}=\{\mathbf{s}|\mathbf{s}\in\mathcal{D},\mathtt{cos}(\mathbf{q},\mathbf{s})\geq\theta\}$. A vector $\mathbf{s}$ is called $\theta$-similar to the query $\mathbf{q}$ if $\mathtt{cos}(\mathbf{q},\mathbf{s})\geq\theta$ and the score of $\mathbf{s}$ is the value $\mathtt{cos}(\mathbf{q},\mathbf{s})$ when $\mathbf{q}$ is understood from the context.

Observe that cosine similarity is insensitive to vector normalization. We will therefore assume without loss of generality that the database as well as query consist of unit vectors (otherwise, all vectors can be normalized in a pre-processing step).

In the literature, cosine threshold querying is a special case of Cosine Similarity Search (CSS) [67, 4, 51], where other aspects like approximate answers, top-k queries and similarity join are considered. Our work considers specifically CSS with exact, threshold and single-vector queries, which is the case of interest to many applications.

Because of the unit-vector assumption, the scoring function $\mathtt{cos}$ computes the dot product $\mathbf{q}\cdot\mathbf{s}$. Without the unit-vector assumption, the problem is equivalent to inner product threshold querying, which is of interest in its own right. Related work on cosine and inner product similarity search is summarized in Section 5.

In this paper we develop novel techniques for the efficient evaluation of cosine threshold queries. We take as a starting point the well-known Threshold Algorithm (TA), by Fagin et al. [29], because of its simplicity, wide applicability, and optimality guarantees. We review the classic $\mathsf{TA}$ in Appendix A.

A TA-like baseline index and algorithm and its shortcomings. The $\mathsf{TA}$ algorithm can be easily adapted to our setting, yielding a first-cut approach to processing cosine threshold queries. We describe how this is done and refer to the resulting index and algorithm as the TA*-like baseline*. Note first that cosine threshold queries use $\mathtt{cos}(\mathbf{q},\mathbf{s})$, which can be viewed as a particular family of functions $F(\mathbf{s})=\mathbf{s}\cdot\mathbf{q}$ parameterized by $\mathbf{q}$, that are monotonic in $\mathbf{s}$ for unit vectors. However, TA produces the vectors with the top-k scores according to $F(\mathbf{s})$, whereas cosine threshold queries return all $\mathbf{s}$ whose score exceeds the threshold $\theta$. We will show how this difference can be overcome straightforwardly.

A baseline index and algorithm inspired by $\mathsf{TA}$ can answer cosine threshold queries exactly without a full scan of the vector database for each query. In addition, the baseline algorithm enjoys the same instance optimality guarantee as the original $\mathsf{TA}$. This baseline is created as follows. First, identically to the TA, the baseline index consists of one sorted list for each of the $d$ dimensions. In particular, the $i$-th sorted list has pairs $(\mathsf{ref}(\mathbf{s}),\mathbf{s}[i])$, where $\mathsf{ref}(\mathbf{s})$ is a reference to the vector $\mathbf{s}$ and $\mathbf{s}[i]$ is its value on the $i$-th dimension. The list is sorted in descending order of $\mathbf{s}[i]$.111There is no need to include pairs with zero values in the list.

Next, the baseline, like the TA, proceeds into a gathering phase during which it collects a complete set of references to candidate result vectors. The TA shows that gathering can be achieved by reading the $d$ sorted lists from top to bottom and terminating early when a stopping condition is finally satisfied. The condition guarantees that any vector that has not been seen yet has no chance of being in the query result. The baseline makes a straightforward change to the TA’s stopping condition to adjust for the difference between the TA’s top-k requirement and the threshold requirement of the cosine threshold queries. In particular, in each round the baseline algorithm has read the first $b$ entries of each index. (Initially it is $b=0$.) If it is the case that $\mathtt{cos}(\mathbf{q},[L_{1}[b],\ldots,$ $L_{d}[b]])<\theta$ then it is guaranteed that the algorithm has already read (the references to) all the possible candidates and thus it is safe to terminate the gathering phase, see Figure 1 for an example. Every vector $\mathbf{s}$ that appears in the $j$-th entry of a list for $j<b$ is a candidate.

In the next phase, called the verification phase, the baseline algorithm (again like TA) retrieves the candidate vectors from the database and checks which ones actually score above the threshold.

For inner product queries, the baseline algorithm’s gathering phase benefits from the same $d\cdot\mathsf{OPT}$ instance optimality guarantee as the TA. Namely, the gathering phase will access at most $d\cdot\mathsf{OPT}$ entries, where $\mathsf{OPT}$ is the optimal index access cost. More specifically, the notion of $\mathsf{OPT}$ is the minimal number of sequential accesses of the sorted inverted index during the gathering phase for any $\mathsf{TA}$-like algorithm applied to the specific query and index instance.

There is an obvious optimization: Only the $k$ dimensions that have non-zero values in the query vector $\mathbf{q}$ should participate in query processing – this leads to a $k\cdot\mathsf{OPT}$ guarantee for inner product queries.222This optimization is equally applicable to the TA’s problem: Scan only the lists that correspond to dimensions that actually affect the function $F$. But even this guarantee loses its practical value when $k$ is a large number. In the mass spectrometry scenario $k$ is $\sim$100. In document similarity and image similarity cases it is even higher.

For cosine threshold queries, the $k\cdot\mathsf{OPT}$ guarantee no longer holds. The baseline fails to utilize the unit vector constraint to reach the stopping condition faster, resulting in an unbounded gap from $\mathsf{OPT}$ because of the unnecessary accesses (see Appendix C).333Notice, the unit vector constraint enables inference about the collective weight of the unseen coordinates of a vector. Furthermore, the baseline fails to utilize the skewing of the values in the vector’s coordinates (both of the database’s vectors and of the query vector) and the linearity of the similarity function. Intuitively, if the query’s weight is concentrated on a few coordinates, the query processing should overweight the respective lists and may, thus, reach the stopping condition much faster than reading all relevant lists in tandem.

We retain the baseline’s index and the gathering-verification structure which characterizes the family of $\mathsf{TA}$-like algorithms. The decision to keep the gathering and verification stages separate is discussed in Section 2. We argue that this algorithmic structure is appropriate for cosine threshold queries, because further optimizations that would require merging the two phases are only likely to yield marginal benefits. Within this framework, we reconsider

Traversal strategy optimization: A traversal strategy determines the order in which the gathering phase proceeds in the lists. In particular, we allow the gathering phase to move deeper in some lists and less deep in others. For example, the gathering phase may have read at some point $b_{1}=106$ entries from the first list, $b_{2}=523$ entries from the second list, etc. Multiple traversal strategies are possible and, generally, each traversal strategy will reach the stopping condition with a different configuration of $[b_{1},b_{2},\ldots,b_{n}]$. The traversal strategy optimization problem asks that we efficiently identify a traversal path that minimizes the access cost $\sum_{i=1}^{d}b_{i}$. To enable such optimization, we will allow lightweight additions to the baseline index. 2. 2.

Stopping condition optimization: We reconsider the stopping condition so that it takes into account (a) the specifics of the $\mathtt{cos}$ function and (b) the unit vector constraint. Moreover, since the stopping condition is tested frequently during the gathering phase, it has to to be evaluated very efficiently. Notice that optimizing the stopping condition is independent of the traversal strategy or skewness assumptions about the data.

**Contributions and summary of results. **

• •

We present a stopping condition for early termination of the index traversal (Section 3). We show that the stopping condition is complete and tight, informally meaning that (1) for any traversal strategy, the gathering phase will produce a candidate set containing all the vectors $\theta$-similar to the query, and (2) the gathering terminates as soon as no more $\theta$-similar vectors can be found (Theorem 3.3). In contrast, the stopping condition of the (TA-inspired) baseline is complete but not tight (Theorem C.1, Appendix C). The proposed stopping condition takes into account that all database vectors are normalized and reduces the problem to solving a special quadratic program (Equation 1) that guarantees both completeness and tightness. While the new stopping condition prunes further the set of candidates, it can also be efficiently computed in $\mathcal{O}(\log d)$ time using incremental maintenance techniques.

• •

We introduce the hull-based traversal strategies that exploit the skewness of the data (Section 4). In particular, skewness implies that each sorted list $L_{i}$ is “mostly convex”, meaning that the shape of $L_{i}$ is approximately the lower convex hull constructed from the set of points of $L_{i}$. This technique is quite general, as it can be extended to the class of decomposable functions which have the form $F(\mathbf{s})=f_{1}(\mathbf{s}[1])+\ldots+f_{d}(\mathbf{s}[d])$ where each $f_{i}$ is non-decreasing.444The inner product threshold problem is the special case where $f_{i}(\mathbf{s}[i])=q_{i}\cdot\mathbf{s}[i]$. Consequently, we provide the following optimality guarantee for inner product threshold queries: The number of accesses executed by the gathering phase (i.e., $\sum_{i=1}^{d}b_{i}$) is at most $\mathsf{OPT}+c$ (Theorem 4.7 and Corollary 4.10), where $\mathsf{OPT}$ is the number of accesses by the optimal strategy and $c$ is the max distance between two vertices in the lower convex hull. Experiments show that in multiple real-world cases, $c$ is a very small fraction of $\mathsf{OPT}$.

• •

Despite the fact that cosine and its tight stopping condition are not decomposable, we show that the hull-based strategy can be adapted to cosine threshold queries by approximating the tight stopping condition with a carefully chosen decomposable function. We show that when the approximation is at most $\epsilon$-away from the actual value, the access cost is at most $\mathsf{OPT}(\theta-\epsilon)+c$ (Theorem 4.12) where $\mathsf{OPT}(\theta-\epsilon)$ is the optimal access cost on the same query $\mathbf{q}$ with the threshold lowered by $\epsilon$ and $c$ is a constant similar to the above decomposable cases. Experiments show that the adjustment $\epsilon$ is very small in practice, e.g., 0.1. We summarize these new results in Table 1.

The paper is organized as follows. We introduce the algorithmic framework and basic definitions in Section 2. Section 3 and 4 discuss the technical developments as we mentioned above. Finally, we discuss related work in Section 5 and conclude in Section 6.

2 Algorithmic Framework

In this section, we present a Gathering-Verification algorithmic framework to facilitate optimizations in different components of an algorithm with a $\mathsf{TA}$-like structure. We start with notations summarized in Table 2.

To support fast query processing, we build an index for the database vectors similar to the original $\mathsf{TA}$. The basic index structure consists of a set of 1-dimensional sorted lists (a.k.a inverted lists in web search [14]) where each list corresponds to a vector dimension and contains vectors having non-zero values on that dimension, as mentioned earlier in Section 1. Formally, for each dimension $i$, $L_{i}$ is a list of pairs $\{(\mathsf{ref}(\mathbf{s}),\mathbf{s}[i])\ |\ \mathbf{s}\in\mathcal{D}\land\mathbf{s}[i]>0\}$ sorted in descending order of $\mathbf{s}[i]$ where $\mathsf{ref}(\mathbf{s})$ is a reference to the vector $\mathbf{s}$ and $\mathbf{s}[i]$ is its value on the $i$-th dimension. In the interest of brevity, we will often write $(\mathbf{s},\mathbf{s}[i])$ instead of $(\mathsf{ref}(\mathbf{s}),\mathbf{s}[i])$. As an example in Figure 1, the list $L_{1}$ is built for the first dimension and it includes 4 entries: $(\mathbf{s}_{1},0.8)$, $(\mathbf{s}_{5},0.7)$, $(\mathbf{s}_{3},0.3)$, $(\mathbf{s}_{4},0.2)$ because $\mathbf{s}_{1}$, $\mathbf{s}_{5}$, $\mathbf{s}_{3}$ and $\mathbf{s}_{4}$ have non-zero values on the first dimension.

Next, we show the Gathering-Verification framework (Algorithm 1) that operates on the index structure. The framework includes two phases: the gathering phase and the verification phase.

Gathering phase (line 1 to line 5). The goal of the gathering phase is to collect a complete set of candidate vectors while minimizing the number of accesses to the sorted lists. The algorithm maintains a position vector $\mathbf{b}=(b_{1},\dots,b_{d})$ where each $b_{i}$ indicates the current position in the inverted list $L_{i}$. Initially, the position vector $\mathbf{b}$ is $(0,\dots,0)$. Then it traverses the lists according to a traversal strategy that determines the list (say $L_{i}$) to be accessed next (line 3). Then it advances the pointer $b_{i}$ by 1 (line 4) and adds the vector $\mathbf{s}$ referenced in the entry $L_{i}[b_{i}]$ to a candidate pool $\mathcal{C}$ (line 5). The traversal strategy is usually stateful, which means that its decision is made based on information that has been observed up to position $\mathbf{b}$ and its past decisions. For example, a strategy may decide that it will make the next 20 moves along dimension 6 and thus it needs state in order to remember that it has already committed to 20 moves on dimension 6.

The gathering phase terminates once a stopping condition is met. Intuitively, based on the information that has been observed in the index, the stopping condition checks if a complete set of candidates has already been found.

Next, we formally define stopping conditions and traversal strategies. As mentioned above, the input of the stopping condition and the traversal strategy is the information that has been observed up to position $\mathbf{b}$, which is formally defined as follows.

Definition 2.1**.**

Let $\mathbf{b}$ be a position vector on the inverted index $\{L_{i}\}_{1\leq i\leq d}$ of a database $\mathcal{D}$. The partial observation at $\mathbf{b}$, denoted as $\mathcal{L}(\mathbf{b})$, is a collection of lists $\{\hat{L}_{i}\}_{1\leq i\leq d}$ where for every $1\leq i\leq d$, $\hat{L_{i}}=[L_{i}[1],\dots,L_{i}[b_{i}]]$.

Definition 2.2**.**

Let $\mathcal{L}(\mathbf{b})$ be a partial observation and $\mathbf{q}$ be a query with similarity threshold $\theta$. A stopping condition is a boolean function $\varphi(\mathcal{L}(\mathbf{b}),\mathbf{q},\theta)$ and a traversal strategy is a function $\mathcal{T}(\mathcal{L}(\mathbf{b}),\mathbf{q},\theta)$ whose domain is $[d]$555$[d]$ is the set $\{1,\dots,d\}$. When clear from the context, we denote them simply by $\varphi(\mathbf{b})$ and $\mathcal{T}(\mathbf{b})$ respectively.

Verification phase (line 6). The verification phase examines each candidate vector $\mathbf{s}$ seen in the gathering phase to verify whether $\mathtt{cos}(\mathbf{q},\mathbf{s})\geq\theta$ by accessing the database. Various techniques [67, 5, 51] have been proposed to speed up this process. Essentially, instead of accessing all the $d$ dimensions of each $\mathbf{s}$ and $\mathbf{q}$ to compute exactly the cosine similarity, these techniques decide $\theta$-similarity by performing a partial scan of each candidate vector. We review these techniques, which we refer to as partial verification, in Appendix B. Additionally, as a novel contribution, we show that in the presence of data skewness, partial verification can have a near-constant performance guarantee (Theorem B.4) for each candidate.

Theorem 2.3**.**

(Informal) For most skewed vectors, $\theta$-similarity can be computed at constant time.

**Remark on optimizing the gathering phase. ** Due to these optimization techniques, the number of sequential accesses performed during the gathering phase becomes the dominating factor of the overall running time. This reason behind is that the number of sequential accesses is strictly greater than the number of candidates that need to be verified so reducing the sequential access cost also results in better performance of the verification phase. In practice, we observed that the sequential cost is indeed dominating: for 1,000 queries on 1.2 billion vectors with similarity threshold 0.6, the sequential gathering time is 16 seconds and the verification time is only 4.6 seconds. Such observation justifies our goal of designing a traversal strategy with near-optimal sequential access cost, as the dominant cost concerns the gathering stage.

Remark on the suitability of $\mathsf{TA}$-like algorithms. One may wonder whether algorithms that start the gathering phase NOT from the top of the inverted lists may outperform the best $\mathsf{TA}$-like algorithm. In particular, it appears tempting to start the gathering phase from the point closest to $q_{i}$ in each inverted list and traverse towards the two ends of each list. Appendix E proves why this idea can lead to poor performance. In particular, Appendix E proves that in a general setting, the computation of a tight and complete stopping condition (formally defined in Definition 3.1 and 3.2) becomes np-hard since it needs to take into account constraints from two pointers (forward and backward) for each inverted list. Furthermore, in many applications, the data skewing leads to small savings from pruning the top area of each list, since the top area is sparsely populated - unlike the densely populated bottom area of each list. Thus it is not justified to use an expensive gathering phase algorithm for small savings.

Section 5.1 reviews additional prior work ideas [67, 68] that avoid traversing some top/bottom regions of the inverted index. Such ideas may provide additional optimizations to $\mathsf{TA}$-like algorithms in variations and/or restrictions of the problem (e.g., a restriction that the threshold is very high) and thus they present future work opportunities in closely related problems.

3 Stopping condition

In this section, we introduce a fine-tuned stopping condition that satisfies the tight and complete requirements to early terminate the index traversal.

First, the stopping condition has to guarantee completeness (Definition 3.1), i.e. when the stopping condition $\varphi$ holds on a position $\mathbf{b}$, the candidate set $\mathcal{C}$ must contain all the true results. Note that since the input of $\varphi$ is the partial observation at $\mathbf{b}$, we must guarantee that for all possible databases $\mathcal{D}$ consistent with the partial observation $\mathcal{L}(\mathbf{b})$, the candidate set $\mathcal{C}$ contains all vectors in $\mathcal{D}$ that are $\theta$-similar to the query $\mathbf{q}$. This is equivalent to require that if a unit vector $\mathbf{s}$ is found below position $\mathbf{b}$ (i.e. $\mathbf{s}$ does not appear above $\mathbf{b}$), then $\mathbf{s}$ is NOT $\theta$-similar to $\mathbf{q}$. We formulate this as follows.

Definition 3.1** (Completeness).**

Given a query $\mathbf{q}$ with threshold $\theta$, a position vector $\mathbf{b}$ on index $\{L_{i}\}_{1\leq i\leq d}$ is complete iff for every unit vector $\mathbf{s}$, $\mathbf{s}<L[\mathbf{b}]$ implies $\mathbf{s}\cdot\mathbf{q}<\theta$. A stopping condition $\varphi(\cdot)$ is complete iff for every $\mathbf{b}$, $\varphi(\mathbf{b})=\mathtt{True}$ implies that $\mathbf{b}$ is complete.

The second requirement of the stopping condition is tightness. It is desirable that the algorithm terminates immediately once the candidate set $\mathcal{C}$ contains a complete set of candidates, such that no additional unnecessary access is made. This can reduce not only the number of index accesses but also the candidate set size, which in turn reduces the verification cost. Formally,

Definition 3.2** (Tightness).**

A stopping condition $\varphi(\cdot)$ is tight iff for every complete position vector $\mathbf{b}$, $\varphi(\mathbf{b})=\mathtt{True}$.

It is desirable that a stopping condition be both complete and tight. However, as we shown in Appendix C, the baseline stopping condition \varphi_{\mathsf{BL}}=\big{(}\mathbf{q}\cdot L[\mathbf{b}]<\theta\big{)} is complete but not tight as it does not capture the unit vector constraint to terminate as soon as no unseen unit vector can satisfy $\mathbf{s}\cdot\mathbf{q}\geq\theta$. Next, we present a new stopping condition that is both complete and tight.

To guarantee tightness, one can check at every snapshot during the traversal whether the current position vector $\mathbf{b}$ is complete and stop once the condition is true. However, directly testing the completeness is impractical since it is equivalent to testing whether there exists a real vector $\mathbf{s}=(s_{1},\dots,s_{d})$ that satisfies the following following set of quadratic constraints:

[TABLE]

We denote by $\mathbf{C}(\mathbf{b})$ (or simply $\mathbf{C}$) the set of $\mathbb{R}^{d}$ points defined by the above constraints. The set $\mathbf{C}(\mathbf{b})$ is infeasible (i.e. there is no satisfying $\mathbf{s}$) if and only if $\mathbf{b}$ is complete, but directly testing the feasibility of $\mathbf{C}(\mathbf{b})$ requires an expensive call to a quadratic programming solver. Depending on the implementation, the running time can be exponential or of high-degree polynomial [13]. We address this challenge by deriving an equivalently strong stopping condition that guarantees tightness and is efficiently testable:

Theorem 3.3**.**

Let $\tau$ be the solution of the equation $\sum_{i=1}^{d}\min\{q_{i}\cdot\tau,L_{i}[b_{i}]\}^{2}=1$ and

[TABLE]

*called the max-similarity. The stopping condition $\varphi_{\mathsf{TC}}(\mathbf{b})=\left(\mathsf{MS}(L[\mathbf{b}])<\theta\right)$ is tight and complete. *

Proof 3.4**.**

The tight and complete stopping condition is obtained by applying the Karush-Kuhn-Tucker (KKT) conditions [45] for solving nonlinear programs. We first formulate the set of constraints in (1) as an optimization problem over $\mathbf{s}$:

[TABLE]

So checking whether $\mathbf{C}$ is feasible is equivalent to verifying whether the maximal $\sum_{i=1}^{d}s_{i}\cdot q_{i}$ is at least $\theta$. So it is sufficient to show that $\sum_{i=1}^{d}s_{i}\cdot q_{i}$ is maximized when $s_{i}=\min\{q_{i}\cdot\tau,L_{i}[b_{i}]\}$ as specified above.

The KKT conditions of the above maximization problem specify a set of necessary conditions that the optimal $\mathbf{s}$ needs to satisfy. More precisely, let

[TABLE]

be the Lagrangian of (3) where $\lambda\in\mathbb{R}$ and $\mathbf{\mu}\in\mathbb{R}^{d}$ are the Lagrange multipliers. Then,

Lemma 3.5** (derived from KKT).**

The optimal $\mathbf{s}$ in (3) satisfies the following conditions:

[TABLE]

in addition to the constraints in (3) (called the Primal feasibility conditions).

By the Complementary slackness condition, for every $i$, if $\mu_{i}\neq 0$ then $s_{i}=L_{i}[b_{i}]$. If $\mu_{i}=0$, then from the Stationarity condition, we know that for every $i$, $q_{i}+\mu_{i}-2\lambda\cdot s_{i}=0$ so $s_{i}=q_{i}/2\lambda$. Thus, the value of $s_{i}$ is either $L_{i}[b_{i}]$ or $q_{i}/2\lambda$.

If $L_{i}[b_{i}]<q_{i}/2\lambda$ then since $s_{i}\leq L_{i}[b_{i}]$, the only possible case is $s_{i}=L_{i}[b_{i}]$. For the remaining dimensions, the objective function $\sum_{i=1}^{d}s_{i}\cdot q_{i}$ is maximized when each $s_{i}$ is proportional to $q_{i}$, so $s_{i}=q_{i}/2\lambda$. Combining these two cases, we have $s_{i}=\min\{q_{i}/2\lambda,L_{i}[b_{i}]\}.$

Thus, for the $\lambda$ that satisfies $\sum_{i=1}^{d}\min\{q_{i}/2\lambda,L_{i}[b_{i}]\}^{2}=1$, the objective function $\sum_{i=1}^{d}s_{i}\cdot q_{i}$ is maximized when $s_{i}=\min\{q_{i}/2\lambda,L_{i}[b_{i}]\}$ for every $i$. The theorem is obtained by letting $\tau=1/2\lambda$.

**Remark of $\varphi_{\mathsf{TC}}$. ** The tight stopping condition $\varphi_{\mathsf{TC}}$ computes the vector $\mathbf{s}$ below $L(\mathbf{b})$ with the maximum cosine similarity $\mathsf{MS}(L[\mathbf{b}])$ with the query $\mathbf{q}$. At the beginning of the gathering phase, $b_{i}=0$ for every $i$ so $\mathsf{MS}(L[\mathbf{b}])=1$ as $\mathbf{s}$ is not constrained. The cosine score is maximized when $\mathbf{s}=\mathbf{q}$ where $\tau=1$. During the gathering phase, as $b_{i}$ increases, the upper bound $L_{i}[b_{i}]$ of each $s_{i}$ decreases. When $L_{i}[b_{i}]<q_{i}$ for some $i$, $s_{i}$ can no longer be $q_{i}$. Instead, $s_{i}$ equals $L_{i}[b_{i}]$, the rest of $\mathbf{s}$ increases proportional to $\mathbf{q}$ and $\tau$ increases. During the traversal, the value of $\tau$ monotonically increases and the score $\mathbf{s}(L[\mathbf{b}])$ monotonically decreases. This is because the space for $\mathbf{s}$ becomes more constrained by $L(\mathbf{b})$ as the pointers move deeper in the inverted lists.

Testing the tight and complete condition $\varphi_{\mathsf{TC}}$ requires solving $\tau$ in Theorem (3.3), for which a direct application of the bisection method takes $\mathcal{O}(d)$ time. We show a novel efficient algorithm (Appendix D) based on incremental maintenance which takes only $\mathcal{O}(\log d)$ time for each test of $\varphi_{\mathsf{TC}}$.

Theorem 3.6**.**

The stopping condition $\varphi_{\mathsf{TC}}(\mathbf{b})$ can be incrementally computed in $\mathcal{O}(\log d)$ time.

4 Near-Optimal Traversal Strategy

Given the inverted lists index and a query, there can be many stopping positions that are both complete and tight. To optimize the performance, we need a traversal strategy that reaches one such position as fast as possible. Specifically, the goal is to design a traversal strategy $\mathcal{T}$ that minimizes $|\mathbf{b}|=\sum_{i=1}^{d}b_{i}$ where $\mathbf{b}$ is the first position vector satisfying the tight and complete stopping condition if $\mathcal{T}$ is followed. Minimizing $|\mathbf{b}|$ also reduces the number of collected candidates, which in turn reduces the cost of the verification phase. We call $|\mathbf{b}|$ the access cost of the strategy $\mathcal{T}$. Formally,

Definition 4.1** (Access Cost).**

Given a traversal strategy $\mathcal{T}$, we denote by $\{\mathbf{b}_{i}\}_{i\geq 0}$ the sequence of position vectors obtained by following $\mathcal{T}$. The access cost of $\mathcal{T}$, denoted by $\mathsf{cost}(\mathcal{T})$, is the minimal $k$ such that $\varphi_{\mathsf{TC}}(\mathbf{b}_{k})=\mathtt{True}$. Note that $\mathsf{cost}(\mathcal{T})$ also equals $|\mathbf{b}_{k}|$.

Definition 4.2** (Instance Optimality).**

Given a database $\mathcal{D}$ with inverted lists $\{L_{i}\}_{1\leq i\leq d}$, a query vector $\mathbf{q}$ and a threshold $\theta$, the optimal access cost $\mathsf{OPT}(\mathcal{D},\mathbf{q},\theta)$ is the minimum $\sum_{i=1}^{d}b_{i}$ for position vectors $\mathbf{b}$ such that $\varphi_{\mathsf{TC}}(\mathbf{b})=\mathtt{True}$. When it is clear from the context, we simply denote $\mathsf{OPT}(\mathcal{D},\mathbf{q},\theta)$ as $\mathsf{OPT}(\theta)$ or $\mathsf{OPT}$.

At a position $\mathbf{b}$, a traversal strategy makes its decision locally based on what has been observed in the inverted lists up to that point, so the capability of making globally optimal decisions is limited. As a result, traversal strategies are often designed as simple heuristics, such as the lockstep strategy in the baseline approach. The lockstep strategy has a $d\cdot\mathsf{OPT}$ near-optimal bound which is loose in the high-dimensionality setting.

In this section, we present a traversal strategy for cosine threshold queries with tighter near-optimal bound by taking into account that the index values are skewed in many realistic scenarios. We approach the (near-)optimal traversal strategy in two steps.

First, we consider the simplified case with the unit-vector constraint ignored so that the problem is reduced to inner product queries. We propose a general traversal strategy that relies on convex hulls pre-computed from the inverted lists during indexing. During the gathering phase, these convex hulls are accessed as auxiliary data during the traversal to provide information on the increase/decrease rate towards the stopping condition. The hull-based traversal strategy not only makes fast decisions (in $\mathcal{O}(\log d)$ time) but is near-optimal (Corollary 4.10) under a reasonable assumption. In particular, we show that if the distance between any two consecutive convex hull vertices of the inverted lists is bounded by a constant $c$, the access cost of the strategy is at most $\mathsf{OPT}+c$. Experiments on real data show that this constant is small in practice.

The hull-based traversal strategy is quite general, as it applies to a large class of functions beyond inner product called the decomposable functions, which have the form $\sum_{i=1}^{d}f_{i}(s_{i})$ where each $f_{i}$ is a non-decreasing real function of a single dimension $s_{i}$. Obviously, for a fixed query $\mathbf{q}$, the inner product $\mathbf{q}\cdot\mathbf{s}$ is a special case of decomposable functions, where each $f_{i}(s_{i})=q_{i}\cdot s_{i}$. We show that the near-optimality result for inner product queries can be generalized to any decomposable function (Theorem 4.7).

Next, in Section 4.4, we consider the cosine queries by taking the normalization constraint into account. Although the function $\mathsf{MS}(\cdot)$ used in the tight stopping condition $\varphi_{\mathsf{TC}}$ is not decomposable so the same technique cannot be directly applied, we show that the hull-based strategy can be adapted by approximating $\mathsf{MS}(\cdot)$ with a decomposable function. In addition, we show that with a properly chosen approximation, the hull-based strategy is near-optimal with a small adjustment to the input threshold $\theta$, meaning that the access cost is bounded by $\mathsf{OPT}(\theta-\epsilon)+c$ for a small $\epsilon$ (Theorem 4.12). Under the same experimental setting, we verify that $\epsilon$ is indeed small in practice.

4.1 Decomposable Functions

We start with defining the decomposable functions for which the hull-based traversal strategies can be applied:

Definition 4.3** (Decomposable Function).**

A decomposable function $F(\mathbf{s})$ is a $d$-dimensional real function where $F(\mathbf{s})=\sum_{i=1}^{d}f_{i}(s_{i})$ and each $f_{i}$ is a non-decreasing real function.

Given a decomposable function $F$, the corresponding stopping condition is called a decomposable condition, which we define next.

Definition 4.4** (Decomposable Condition).**

A decomposable condition $\varphi_{F}$ is a boolean function \varphi_{F}(\mathbf{b})=\big{(}F(L[\mathbf{b}])<\theta\big{)} where $F$ is a decomposable function and $\theta$ is a fixed threshold.

When the unit vector constraint is lifted, the decomposable condition is tight and complete for any scoring function $F$ and threshold $\theta$. As a result, the goal of designing a traversal strategy for $F$ is to have the access cost as close as possible to $\mathsf{OPT}$ when the stopping condition is $\varphi_{F}$.

4.2 The max-reduction traversal strategy

To illustrate the high-level idea of the hull-based approach, we start with a simple greedy traversal strategy called the Max-Reduction traversal strategy $\mathcal{T}_{\mathsf{MR}}(\cdot)$. The strategy works as follows: at each snapshot, move the pointer $b_{i}$ on the inverted list $L_{i}$ that results in the maximal reduction on the score $F(L[\mathbf{b}])$. Formally, we define

[TABLE]

where $\mathds{1}_{i}$ is the vector with 1 at dimension $i$ and 0’s else where. Such a strategy is reasonable since one would like $F(L[\mathbf{b}])$ to drop as fast as possible, so that once it is below $\theta$, the stopping condition $\varphi_{F}$ will be triggered and terminate the traversal.

It is obvious that there are instances where the max-reduction strategy can be far from optimal, but is it possible that it is optimal under some assumption? The answer is positive: if for every list $L_{i}$, the values of $f_{i}(L_{i}[b_{i}])$ are decreasing at decelerating rate, then we can prove that its access cost is optimal. We state this ideal assumption next.

Assumption 1** (Ideal Convexity).**

For every inverted list $L_{i}$, let $\Delta_{i}[j]=f_{i}(L_{i}[j])-f_{i}(L_{i}[j+1])$ for $0\leq j<|L_{i}|$.666Recall that $L_{i}[0]=1$. The list $L_{i}$ is ideally convex if the sequence $\Delta_{i}$ is non-increasing, i.e., $\Delta_{i}[j+1]\leq\Delta_{i}[j]$ for every $j$. Equivalently, the piecewise linear function passing through the points $\{(j,f_{i}(L_{i}[j]))\}_{0\leq j\leq|L_{i}|}$ is convex for each $i$. A database $\mathcal{D}$ is ideally convex if every $L_{i}$ is ideally convex.

An example of an inverted list satisfying the above assumption is shown in Figure 2(a). The max-reduction strategy $\mathcal{T}_{\mathsf{MR}}$ is optimal under the ideal convexity assumption:

Theorem 4.5** (Ideal Optimality).**

Given a decomposable function $F$, for every ideally convex database $\mathcal{D}$ and every threshold $\theta$, the access cost of $\mathcal{T}_{\mathsf{MR}}$ is exactly $\mathsf{OPT}$.

We prove Theorem 4.5 with a simple greedy argument (Appendix F): each move of $\mathcal{T}_{\mathsf{MR}}$ always results in the globally maximal reduction in the scoring function as guaranteed by the convexity condition.

4.3 The hull-based traversal strategy

Theorem 4.5 provides a strong performance guarantee but the ideal convexity assumption is usually not true on real datasets. Without the ideal convexity assumption, the strategy suffers from the drawback of making locally optimal but globally suboptimal decisions. The pointer $b_{i}$ to an inverted list $L_{i}$ might never be moved if choosing the current $b_{i}$ only results in a small decrease in the score $F(L[\mathbf{b}])$, but there is a much larger decrease several steps ahead. As a result, the $\mathcal{T}_{\mathsf{MR}}$ strategy has no performance guarantee in general.

In most practical scenarios that we have seen, we can bring the traversal strategy $\mathcal{T}_{\mathsf{MR}}$ to practicality by considering a relaxed version of Assumption 1. Informally, instead of assuming that each list $f_{i}(L_{i})$ forms a convex piecewise linear function, we assume that $f_{i}(L_{i})$ is “mostly” convex, meaning that if we compute the lower convex hull [22] of $f_{i}(L_{i})$, the gap between any two consecutive vertices on the convex hull is small.777We denote by $f_{i}(L_{i})$ the list $[f_{i}(L_{i}[0]),f_{i}(L_{i}[1]),\dots]$ for every $L_{i}$. Intuitively, the relaxed assumption implies that the values at each list are decreasing at “approximately” decelerating speed. It allows list segments that do not follow the overall deceleration trend, as long as their lengths are bounded by a constant. We verified this property in the mass spectrometry dataset as illustrated in Figure 3, a document dataset, and an image dataset (Figure 11 and 12 in Appendix I).

Assumption 2** (Near-Convexity).**

For every inverted list $L_{i}$, let $\mathsf{H}_{i}$ be the lower convex hull of the set of 2-D points $\{(j,f_{i}(L_{i}[j]))\}_{0\leq j\leq|L_{i}|}$ represented by a set of indices $\mathsf{H}_{i}=\{j_{1},\dots,j_{n}\}$ where for each $1\leq k\leq n$, $(j_{k},f_{i}(L_{i}[j_{k}]))$ is a vertex of the convex hull. The list $L_{i}$ is near-convex if for every $k$, $j_{k+1}-j_{k}$ is upper-bounded by some constant $c$. A database $\mathcal{D}$ is near-convex if every inverted list $L_{i}$ is near-convex with the same constant $c$, which we refer to as the convexity constant.

Example 4.6**.**

*Intuitively, the near-convexity assumption captures the case where each $f_{i}(L_{i})$ is decreasing with approximately decelerating speed, so the number of points between two convex hull vertices should be small. For example, when $f_{i}$ is a linear function, the list $L_{i}$ shown in Figure 2(b) is near-convex with convexity constant 2 since there is at most 1 point between each pair of consecutive vertices of the convex hull (dotted line). In the ideal case shown in Figure 2(a), the constant is 1 when the decrease between successive values is strictly decelerating. *

Imitating the max-reduction strategy, for every pair of consecutive indices $j_{k},j_{k+1}$ in $\mathsf{H}_{i}$ and for every index $j\in[j_{k},j_{k+1})$, let $\tilde{\Delta}_{i}[j]=\dfrac{f_{i}(L_{i}[j_{k}])-f_{i}(L_{i}[j_{k+1}])}{j_{k+1}-j_{k}}$. Since the $(j_{k},f_{i}(L_{i}[j_{k}]))$’s are vertices of a lower convex hull, each sequence $\tilde{\Delta}_{i}$ is non-decreasing. Then the hull-based traversal strategy is simply defined as

[TABLE]

Remark on data structures. In a practical implementation, to answer queries with scoring function $F$ using the hull-based strategy, the lower convex hulls need to be ready before the traversal starts. If $F$ is a general function unknown a priori, the convex hulls need to be computed online which is not practical. Fortunately, when $F$ is the inner product $F(\mathbf{s})=\mathbf{q}\cdot\mathbf{s}$ parameterized by the query $\mathbf{q}$, each convex hull $\mathsf{H}_{i}$ is exactly the convex hull for the points $\{(j,L_{i}[j])\}_{0\leq i\leq|L_{i}|}$ from $L_{i}$. This is because the slope from any two points $(j,f_{i}(L_{i}[j]))$ and $(k,f_{i}(L_{i}[k]))$ is $\dfrac{q_{i}L_{i}[j]-q_{i}L_{i}[k]}{j-k}$, which is exactly the slope from $(j,L_{i}[j])$ and $(k,L_{i}[k])$ multiplied by $q_{i}$. So by using the standard convex hull algorithm [22], $\mathsf{H}_{i}$ can be pre-computed in $\mathcal{O}(|L_{i}|)$ time. Then the set of the convex hull vertices $\mathsf{H}_{i}$ can be stored as inverted lists and accessed for computing the $\tilde{\Delta}_{i}$’s during query processing. In the ideal case, $\mathsf{H}_{i}$ can be as large as $|L_{i}|$ but is much smaller in practice.

Moreover, during the traversal using the strategy $\mathcal{T}_{\mathsf{HL}}$, choosing the maximum $\tilde{\Delta}_{i}[b_{i}]$ at each step can be done in $\mathcal{O}(\log d)$ time using a max heap. This satisfies the requirement that the traversal strategy is efficiently computable.

**Near-optimality results. **We show that the hull-based strategy $\mathcal{T}_{\mathsf{HL}}$ is near-optimal under the near-convexity assumption.

Theorem 4.7**.**

Given a decomposable function $F$, for every near-convex database $\mathcal{D}$ and every threshold $\theta$, the access cost of $\mathcal{T}_{\mathsf{HL}}$ is strictly less than $\mathsf{OPT}+c$ where $c$ is the convexity constant.

When the assumption holds with a small convexity constant, this near-optimality result provides a much tighter bound compared to the $d\cdot\mathsf{OPT}$ bound in the $\mathsf{TA}$-inspired baseline. This is achieved under data assumption and by keeping the convex hulls as auxiliary data structure, so it does not contradict the lower bound results on the approximation ratio [29].

Proof 4.8**.**

Let $\mathcal{B}=\{\mathbf{b}_{i}\}_{i\geq 0}$ be the sequence of position vectors generated by $\mathcal{T}_{\mathsf{HL}}$. We call a position vector $\mathbf{b}$ a boundary position if every $b_{i}$ is the index of a vertex of the convex hull $\mathsf{H}_{i}$. Namely, $b_{i}\in\mathsf{H}_{i}$ for every $i\in[d]$. Notice that if we break ties consistently during the traversal of $\mathcal{T}_{\mathsf{HL}}$, then in between every pair of consecutive boundary positions $\mathbf{b}$ and $\mathbf{b}^{\prime}$ in $\mathcal{B}$, $\mathcal{T}_{\mathsf{HL}}(\mathbf{b})$ will always be the same index. We call the subsequence positions $\{\mathbf{b}_{i}\}_{l\leq i<r}$ of $\mathcal{B}$ where $\mathbf{b}_{l}=\mathbf{b}$ and $\mathbf{b}_{r}=\mathbf{b}^{\prime}$ a segment with boundaries $(\mathbf{b}_{l},\mathbf{b}_{r})$. We show the following lemma.

Lemma 4.9**.**

For every boundary position vector $\mathbf{b}$ generated by $\mathcal{T}_{\mathsf{HL}}$, we have $F(L[\mathbf{b}])\leq F(L[\mathbf{b}^{*}])$ for every position vector $\mathbf{b}^{*}$ where $|\mathbf{b}^{*}|=|\mathbf{b}|$.

Intuitively, the above lemma says that if the traversal of $\mathcal{T}_{\mathsf{HL}}$ reaches a boundary position $\mathbf{b}$, then the score $F(L[\mathbf{b}])$ is the minimal possible score obtained by any traversal sequence of at most $|\mathbf{b}|$ steps. We prove Lemma 4.9 by generalizing the greedy argument in the proof of Theorem 4.5. More details can be found in Appendix G.

Lemma 4.9 is sufficient for Theorem 4.7 because of the following. Suppose $\mathbf{b}_{\mathsf{stop}}$ is the stopping position in $\mathcal{B}$, which means that $\mathbf{b}_{\mathsf{stop}}$ is the first position in $\mathcal{B}$ that satisfies $\varphi_{F}$ and the access cost is $|\mathbf{b}_{\mathsf{stop}}|$. Let $\{\mathbf{b}_{i}\}_{l\leq i<r}$ be the segment that contains $\mathbf{b}_{\mathsf{stop}}$. Given Lemma 4.9, Theorem 4.7 holds trivially if $\mathbf{b}_{\mathsf{stop}}=\mathbf{b}_{l}$. It remains to consider the case $\mathbf{b}_{\mathsf{stop}}\neq\mathbf{b}_{l}$. Since the traversal does not stop at $\mathbf{b}_{l}$, we have $F(L[\mathbf{b}_{l}])\geq\theta$. By Lemma 4.9, $\mathbf{b}_{l}$ is the position with minimal $F(L[\cdot])$ obtained in $|\mathbf{b}_{l}|$ steps so $|\mathbf{b}_{l}|\leq\mathsf{OPT}$. Since $|\mathbf{b}_{\mathsf{stop}}|-|\mathbf{b}_{l}|<|\mathbf{b}_{r}|-|\mathbf{b}_{l}|\leq c$, we have that $|\mathbf{b}_{\mathsf{stop}}|<\mathsf{OPT}+c$. We illustrate this in Figure 4.

Since the baseline stopping condition $\varphi_{\mathsf{BL}}$ is tight and complete for inner product queries, one immediate implication of Theorem 4.7 is that

Corollary 4.10**.**

(Informal) The hull-based strategy $\mathcal{T}_{\mathsf{HL}}$ for inner product queries is near-optimal.

Verifying the assumption. We demonstrate the practical impact of the near-optimality result in real mass spectrometry datasets. The same experiment is repeated on a document and an image dataset in Appendix I. The near-convexity assumption requires that the gap between any two consecutive convex hull vertices has bounded size, which is hard to achieve in general. According to the proof of Theorem 4.7, for a given query, the difference from the optimal access cost is at most the size of the gap between the two consecutive convex hull vertices containing the last move of the strategy (the $\mathbf{b}_{l}$ and $\mathbf{b}_{r}$ in Figure 4). The size of this gap can be much smaller than the global convexity constant $c$, so the overall precision can be much better in practice. We verify this by running a set of 1,000 real queries on the dataset888https://proteomics2.ucsd.edu/ProteoSAFe/index.jsp.

The gap size is 163.04 in average, which takes only 1.3% of the overall access cost of traversing the indices. This indicates that the near-optimality guarantee holds in the mass spectrometry dataset. Similar results are obtained in a document and an image dataset, where the gap size takes only 7.9% and 0.4% of the overall access cost respectively.

4.4 The traversal strategy for cosine

Next, we consider traversal strategies which take into account the unit vector constraint posed by the cosine function, which means that the tight and complete stopping condition is $\varphi_{\mathsf{TC}}$ introduced in Section 3. However, since the scoring function $\mathsf{MS}$ in $\varphi_{\mathsf{TC}}$ is not decomposable, the hull-based technique cannot be directly applied. We adapt the technique by approximating the original $\mathsf{MS}$ with a decomposable function $\tilde{F}$. Without changing the stopping condition $\varphi_{\mathsf{TC}}$, the hull-based strategy can then be applied with the convex hull indices constructed with the approximation $\tilde{F}$. In the rest of this section, we first generalize the result in Theorem 4.7 to scoring functions having decomposable approximations and show how the hull-based traversal strategy can be adapted. Next, we show a natural choice of the approximation for $\mathsf{MS}$ with practically tight near-optimal bounds. Finally, we discuss data structures to support fast query processing using the traversal strategy.

We start with some additional definitions.

Definition 4.11**.**

A $d$-dimensional function $F$ is decomposably approximable if there exists a decomposable function $\tilde{F}$, called the decomposable approximation of $F$, and two non-negative constants $\epsilon_{1}$ and $\epsilon_{2}$ such that $\tilde{F}(\mathbf{s})-F(\mathbf{s})\in[-\epsilon_{1},\epsilon_{2}]$ for every vector $\mathbf{s}$.

When applied to a decomposably approximable function $F$, the hull-based traversal strategy $\mathcal{T}_{\mathsf{HL}}$ is adapted by constructing the convex hull indices and the $\{\tilde{\Delta}_{i}\}_{1\leq i\leq d}$ using the approximation $\tilde{F}$. The following can be obtained by generalizing Theorem 4.7:

Theorem 4.12**.**

Given a function $F$ approximable by a decomposable function $\tilde{F}$ with constants $(\epsilon_{1},\epsilon_{2})$, for every near-convex database $\mathcal{D}$ wrt $\tilde{F}$ and every threshold $\theta$, the access cost of $\mathcal{T}_{\mathsf{HL}}$ is strictly less than $\mathsf{OPT}(\theta-\epsilon_{1}-\epsilon_{2})+c$ where $c$ is the convexity constant.

Proof 4.13**.**

Recall that $\mathbf{b}_{l}$ is the last boundary position generated by $\mathcal{T}_{\mathsf{HL}}$ that does not satisfy the tight stopping condition for $F$ (which is $\varphi_{\mathsf{TC}}$ when $F$ is $\mathsf{MS}$) so $F(L[\mathbf{b}_{l}])\geq\theta$. It is sufficient to show that for every vector $\mathbf{b}^{*}$ where $|\mathbf{b}^{*}|=|\mathbf{b}_{l}|$, $F(L[\mathbf{b}^{*}])\geq\theta-\epsilon_{1}-\epsilon_{2}$ so no traversal can stop within $|\mathbf{b}_{l}|$ steps, implying that the final access cost is no more than $|\mathbf{b}_{l}|+c$ which is bounded by $\mathsf{OPT}(\theta-\epsilon_{1}-\epsilon_{2})+c$.

By Lemma 4.9, we know that for every such $\mathbf{b}^{*}$, $\tilde{F}(L[\mathbf{b}^{*}])\geq\tilde{F}(L[\mathbf{b}_{l}])$. By definition of the approximation $\tilde{F}$, we know that $F(L[\mathbf{b}^{*}])\geq\tilde{F}(L[\mathbf{b}^{*}])-\epsilon_{1}$ and $\tilde{F}(L[\mathbf{b}_{l}])\geq F(L[\mathbf{b}_{l}])-\epsilon_{2}$. Combined together, for every $\mathbf{b}^{*}$ where $|\mathbf{b}^{*}|=|\mathbf{b}_{l}|$, we have

[TABLE]

This completes the proof of Theorem 4.12.

**Choosing the decomposable approximation. ** By Theorem 4.12, it is important to choose an approximation $\tilde{F}$ of $\mathsf{MS}$ with small $\epsilon_{1}$ and $\epsilon_{2}$ for a tight near-optimality result. By inspecting the formula (2) of $\mathsf{MS}$, one reasonable choice of $\tilde{F}$ can be obtained by replacing the term $\tau$ with a fixed constant $\tilde{\tau}$. Formally, let

[TABLE]

be the decomposable approximation of $\mathsf{MS}$ where each component is a non-decreasing function $f_{i}(x)=\min\{q_{i}\cdot\tilde{\tau},x\}\cdot q_{i}$ for $i\in[d]$.

Ideally, the approximation is tight if the constant $\tilde{\tau}$ is close to the final value of $\tau$ which is unknown in advance. We argue that when $\tilde{\tau}$ is properly chosen, the approximation parameter $\epsilon_{1}+\epsilon_{2}$ is very small. With a detailed analysis in Appendix H, we obtain the following upper bound of $\epsilon$:

[TABLE]

**Verifying the near-optimality. ** Next, we verify that the above upper bound of $\epsilon$ is small in practice. We ran the same set of queries as in Section 4.3 and show the distribution of $\epsilon$’s upper bounds in Figure 6. We set $\tilde{\tau}=1/\theta$ for all queries so the first term of (6) becomes zero. Note that more aggressive pruning can yield better $\epsilon$, but it is not done here for simplicity. Overall, the fraction of queries with an upper bound $<$0.12 (the sum of the first 3 bars for all $\theta$) is 82.5% and the fraction of queries with $\epsilon>0.16$ is $0.5\%$. Similar to the case with inner product queries, the average of the convexity constant $c$ is 193.39, which is only 4.8% of the overall access cost.

**Remark on data structures. ** Similar to the inner product case, it is necessary that the convex hulls for $\mathcal{T}_{\mathsf{HL}}$ can be efficiently obtained without a full computation when a query comes in. For every $i\in[d]$, we let $\tilde{\mathsf{H}}_{i}$ be the convex hull for the $i$-th component $f_{i}$ of $\tilde{F}$ and $\mathsf{H}_{i}$ be the convex hull constructed directly from the original inverted list $L_{i}$. Next, we show that each $\tilde{\mathsf{H}}_{i}$ can be efficiently obtained from $\mathsf{H}_{i}$ during query time so we only need to pre-compute the $\mathsf{H}_{i}$’s.

We observe that when $L_{i}[b_{i}]\geq q_{i}\cdot\tilde{\tau}$, $f_{i}(L_{i}[b_{i}])$ equals a fixed value $q_{i}^{2}\cdot\tilde{\tau}$ otherwise is proportional to $L_{i}[b_{i}]$. As illustrated in Figure 6 (left), the list of values $\{f_{i}(L_{i}[j])\}_{j\geq 0}$ is essentially obtained by replacing the $L_{i}[j]$’s greater than $q_{i}\cdot\tilde{\tau}$ with $q_{i}\cdot\tilde{\tau}$.

The following can be shown using properties of convex hulls:

Lemma 4.14**.**

For every $i\in[d]$, the convex hull $\tilde{\mathsf{H}}_{i}$ is a subset of $\mathsf{H}_{i}$ where an index $j_{k}$ of $\mathsf{H}_{i}$ is in $\tilde{\mathsf{H}}_{i}$ iff $k=1$ or

[TABLE]

Lemma 4.14 provides an efficient way to obtain each convex hull $\tilde{\mathsf{H}}_{i}$ from the pre-computed $\mathsf{H}_{i}$’s. When a query $\mathbf{q}$ is given, we perform a binary search on each $\mathsf{H}_{i}$ to find the first $j_{k}\in\mathsf{H}_{i}$ that satisfies (7). Then $\tilde{\mathsf{H}}_{i}$ is the set of indices $\{0,j_{k},j_{k+1}\dots\}$. We illustrate the construction in Figure 6 (right).

Suppose that the maximum size of all $\mathsf{H}_{i}$ is $h$. The computation of the $\tilde{\mathsf{H}}_{i}$’s adds an extra $\mathcal{O}(d\log h)$ of overhead to the query processing time, which is insignificant in practice since $h$ is likely to be much smaller than the size of the database.

5 Related work

In this section, we present the main related work and defer the additional related work (e.g., dimensionality reduction, mass spectrometry search, inner product queries) to the appendix (Section K).

5.1 Cosine similarity search

The cosine threshold querying studied in this work is a special case of the cosine similarity search (CSS) problem [10, 4, 5] mentioned in Section 1. We first survey the techniques developed for CSS.

LSH. A widely used technique for cosine similarity search is locality-sensitive hash (LSH) [59, 6, 36, 38, 66]. The main idea of LSH is to partition the whole database into buckets using a series of hash functions such that similar vectors have high probability to be in the same bucket. However, LSH is designed for approximate query processing, meaning that it is not guaranteed to return all the true results. In contrast, this work focuses on exact query processing which returns all the results.

TA-family algorithms. Another technique for cosine similarity search is the family of $\mathsf{TA}$-like algorithms. Those algorithms were originally designed for processing top-k ranking queries that find the top $k$ objects ranked according to an aggregation function (see [37] for a survey). We have summarized the classic $\mathsf{TA}$ algorithm [29], presented a baseline algorithm inspired by it, and explained its shortcomings in Section 1. The Gathering-Verification framework introduced in Section 2 captures the typical structure of the $\mathsf{TA}$-family when applied to our setting.

The variants of TA (e.g., [33, 9, 23, 15]) can have poor or no performance guarantee for cosine threshold queries since they do not fully leverage the data skewness and the unit vector condition. For example, Güntzer et al. developed Quick-Combine [33]. Instead of accessing all the lists in a lockstep strategy, it relies on a heuristic traversal strategy to access the list with the highest rate of changes to the ranking function in a fixed number of steps ahead. It was shown in [30] that the algorithm is not instance optimal. Although the hull-based traversal strategy proposed in this paper roughly follows the same idea, the number of steps to look ahead is variable and determined by the next convex hull vertex. Thus, for decomposable functions, the hull-based strategy makes globally optimal decisions and is near-optimal under the near-convexity assumption, while Quick-Combine has no performance guarantee because of the fixed step size even when the data is near-convex. Other TA variants are discussed in Appendix K.2.

COORD. Teflioudi et al. proposed the COORD algorithm based on inverted lists for CSS [68, 67]. The main idea is to scan the whole lists but with an optimization to prune irrelevant entries using upper/lower bounds of the cosine similarity with the query. Thus, instead of traversing the whole lists starting from the top, it scans only those entries within a feasible range. We can also apply such a pruning strategy to the Gathering-Verification framework by starting the gathering phase at the top of the feasible range. However, there is no optimality guarantee of the algorithm. Also the optimization only works for high thresholds (e.g., 0.95), which are not always the requirement. For example, a common and well-accepted threshold in mass spectrometry search is 0.6, which is a medium-sized threshold, making the effect of the pruning negligible.

Partial verification. Anastasiu and Karypis proposed a technique for fast verification of $\theta$-similarity between two vectors [4] without a full scan of the two vectors. We can apply the same optimization to the verification phase of the Gathering-Verification framework. Additionally, we prove that it has a novel near-constant performance guarantee in the presence of data skewness. See Appendix B.

Other variants. There are several studies focusing on cosine similarity join to find out all pairs of vectors from the database such that their similarity exceeds a given threshold [10, 4, 5]. However, this work is different since the focus is comparing to a given query vector q rather than join. As a result, the techniques in [10, 4, 5] are not directly applicable: (1) The inverted index is built online instead of offline, meaning that at least one full scan of the whole data is required, which is inefficient for search. (2) The index in  [10, 4, 5] is built for a fixed query threshold, meaning that the index cannot be used for answering arbitrary query thresholds as concerned in this work. The theoretical aspects of similarity join were discussed recently in [2, 36].

5.2 Euclidean distance threshold queries

The cosine threshold queries can also be answered by techniques for distance threshold queries (the threshold variant of nearest neighbor search) in Euclidean space. This is because there is a one-to-one mapping between the cosine similarity $\theta$ and the Euclidean distance $r$ for unit vectors, i.e., $r=2\sin(\arccos(\theta)/2)$. Thus, finding vectors that are $\theta$-similar to a query vector is equivalent to finding the vectors whose Euclidean distance is within $r$. Next, we review exact approaches for distance queries while leaving the discussion of approximate approaches in the appendix (Section K.1). There are four main types of techniques for exact approaches: tree-based indexing, pivot-based indexing, clustering, and dimensionality reduction (See Appendix K.5).

Tree-based indexing. Several tree-based indexing techniques (such as R-tree, KD-tree, Cover-tree [11]) were developed for range queries (so they can also be applied to distance queries), see [12] for a survey. However, they are not scalable to high dimensions (say thousands of dimensions as studied in this work) due to the well known dimensionality curse issue [74].

Pivot-based indexing. The main idea is to pre-compute the distances between data vectors and a set of selected pivot vectors. Then during query processing, use triangle inequalities to prune irrelevant vectors [17, 35]. However, it does not scale in high-dimensional space as shown in [17] since it requires a large space to store the pre-computed distances.

Clustering-based (or partitioning-based) methods. The main idea of clustering is to partition the database vectors into smaller clusters of vectors during indexing. Then during query processing, irrelevant clusters are pruned via the triangle inequality [62, 61]. Clustering is an optimization orthogonal to the proposed techniques, as they can be used to process vectors within each cluster to speed up the overall performance.

6 Conclusion

In this work, we proposed optimizations to the index-based, TA-like algorithms for answering the cosine threshold queries, which lie at the core of numerous applications. The novel techniques include a complete and tight stopping condition computable incrementally in $\mathcal{O}(\log d)$ time and a family of convex hull-based traversal strategies with near-optimality guarantees for a larger class of decomposable functions beyond cosine. With these techniques, we show near-optimality first for inner-product threshold queries, then extend the result to the full cosine threshold queries using approximation. These results are significant improvements over a baseline approach inspired by the classic $\mathsf{TA}$ algorithm. In addition, we have verified with experiments on real data the assumptions required by the near-optimality results.

Appendix

Appendix A Review of $\mathsf{TA}$

On a database $\mathcal{D}$ of $d$-dimensional vectors $\{\mathbf{s}_{1},\dots,\mathbf{s}_{n}\}$, given a query monotonic scoring function $F:\mathbb{R}^{d}\mapsto\mathbb{R}$ and a query parameter $k$, the Threshold Algorithm ($\mathsf{TA}$) computes the $k$ database vectors with the highest score $F(\mathbf{s})$. First, $\mathsf{TA}$ preprocesses the vector database by building an inverted index $\{L_{i}\}_{1\leq i\leq d}$ where each $L_{i}$ is an inverted list contains pairs of $(\mathsf{ref}(\mathbf{s}),\mathbf{s}[i])$ where $\mathsf{ref}(\mathbf{s})$ is a reference of the vector $\mathbf{s}$ and $\mathbf{s}[i]$ is the $i$-th dimension $\mathbf{s}$. Each $L_{i}$ is sorted in descending order of $\mathbf{s}[i]$. When a query $(F,k)$ arrives, $\mathsf{TA}$ proceeds as follows. It maintains a pointer $b$ starting from 0 to all the inverted lists and increments $b$ iteratively. At each iteration:

• •

Collect the set $\mathcal{C}$ of candidates of all references in $L_{i}$ upto position $b$ for all dimension $i$. Namely, $\mathcal{C}=\bigcup_{i=1}^{d}\{\mathsf{ref}(\mathbf{s})|(\mathsf{ref}(\mathbf{s}),\mathbf{s}[i])\in L_{i}[0,\dots,b]\}$.

• •

Compute $F_{k}$ the $k$-th highest score $F(\mathbf{s})$ for all $\mathsf{ref}(\mathbf{s})\in\mathcal{C}$ by accessing $\mathbf{s}$ in the database with the reference.

• •

If the score $F_{k}$ is no less than $F(L_{1}[b],\dots,L_{d}[b])$, return the top-$k$ highest score vectors in $\mathcal{C}$; otherwise continue to the next iteration with $b\leftarrow b+1$.

By monotonicity of the function $F$, once the stopping condition is satisfied, it is guaranteed that no vector $\mathbf{s}$ below the pointer $b$ can have $F(\mathbf{s})$ above the current $k$-th highest score. Thus the candidate set $\mathcal{C}$ contains the complete set of all the $k$ highest score vectors in the database.

A.1 The $\mathsf{TA}$-inspired Baseline

Next, we show the baseline stopping condition and traversal strategy inspired by the $\mathsf{TA}$ for cosine threshold queries. As already reviewed in Section 1, the gathering phase stops when the cosine function score at the current position $\mathbf{b}$ is below the input threshold $\theta$. Formally, the baseline stopping condition $\varphi_{\mathsf{BL}}$ is

[TABLE]

To determine the order of accessing the inverted lists, $\mathsf{TA}$ advances all the pointers at equal rate. An obvious optimization is to move only the pointers whose dimension has non-zero values in $\mathbf{q}$. Let $\mathsf{nz}=\{i|i\in[d]\land q_{i}>0\}$ be the list of all non-zero dimensions and let $m=|\mathsf{nz}|$. The baseline traversal strategy $\mathcal{T}_{\mathsf{BL}}$ is the following:

[TABLE]

so that each inverted list of a non-zero dimension is advanced at equal speed. By the classic result of [29],

Theorem A.1**.**

For inner product threshold queries (without the normalization constraint), the access cost of the $\mathsf{TA}$-inspired baseline is at most $m\cdot\mathsf{OPT}$ where $m$ is the query’s number of non-zero dimensions and $\mathsf{OPT}$ is the optimal access cost.

For cosine threshold queries, since the baseline condition $\varphi_{\mathsf{BL}}$ is not tight (Theorem C.1), the same instance optimality bound does not hold.

Appendix B The Verification Phase

Next, we discuss optimizations in the verification phase where each gathered candidate is tested for $\theta$-similarity with the query. The naive approach of verification is to fully access all the non-zero entries of each candidate s to compute the exact similarity score $\mathtt{cos}(\mathbf{q},\mathbf{s})$ and compare it against $\theta$, which takes $\mathcal{O}(d)$ time. Various techniques have been proposed [67, 5, 51] to decide $\theta$-similarity by leveraging only partial information about the candidate vector so the exact computation of $\mathtt{cos}(\mathbf{q},\mathbf{s})$ can be avoided. In this section, we revisit these existing techniques which we call partial verification. In addition, as a novel contribution, we show that in the presence of data skewness, partial verification can have a near-constant performance guarantee (Theorem B.4).

Informally, while a vector $\mathbf{s}$ is scanned, based on what has been observed in $\mathbf{s}$ so far, it is possible to infer that

• (1)

the similarity score $\mathtt{cos}(\textbf{q},\mathbf{s})$ is certainly at least $\theta$ or

• (2)

the similarity score is certainly below $\theta$.

In either case, we can stop without scanning the rest of $\mathbf{s}$ and return an accurate verification result. The problem is formally defined as follows:

Problem B.1** (Partial Verification).**

A partially observed vector $\tilde{\mathbf{s}}$ is a $d$-dimensional vector in $(\mathbb{R}_{\geq 0}\cup\{\bot\})^{d}$. Given a query $\mathbf{q}$ and a partially observed vector $\tilde{\mathbf{s}}$, compute whether for every vector $\mathbf{s}$ where $\mathbf{s}[i]=\tilde{\mathbf{s}}[i]$ for every $\tilde{\mathbf{s}}[i]\neq\bot$, it is $\mathtt{cos}(\mathbf{s},\mathbf{q})\geq\theta$.

Intuitively, a partially observed vector $\tilde{\mathbf{s}}$ contains the entries of a candidate $\mathbf{s}$ either already observed during the gathering phase, or accessed for the actual values during the verification phase. The unobserved dimensions are replaced with a null value $\bot$. We say that a vector $\mathbf{s}$ is compatible with a partially observed one $\tilde{\mathbf{s}}$ if $\mathbf{s}[i]=\tilde{\mathbf{s}}[i]$ for every dimension $i$ where $\tilde{\mathbf{s}}[i]\neq\bot$.

The partial verification problem can be solved by computing an upper and a lower bound of the cosine similarity between $\mathbf{s}$ and $\mathbf{q}$ when $\tilde{\mathbf{s}}$ is observed. We denote by $ub(\tilde{\mathbf{s}})$ and $lb(\tilde{\mathbf{s}})$ the upper bound and lower bound, so $ub(\tilde{\mathbf{s}})=\max_{\mathbf{s}}\{\mathbf{s}\cdot\mathbf{q}\}$ and $lb(\tilde{\mathbf{s}})=\min_{\mathbf{s}}\{\mathbf{s}\cdot\mathbf{q}\}$ where the maximum/minimum are taken over all $\mathbf{s}$ compatible with $\tilde{\mathbf{s}}$. By [67, 5, 51], the upper/lower bounds can be computed as follows:

Lemma B.2**.**

Given a partially observed vector $\tilde{\mathbf{s}}$ and a query vector $\mathbf{q}$,

[TABLE]

Example B.3**.**

Figure 7 shows an example of computing the lower/upper bounds of a partially scanned vector $\mathbf{s}$ with a query $\mathbf{q}$. Assume the first three dimensions have been scanned, then the lower bound and upper bound are computed as follows:

[TABLE]

If the threshold $\theta$ is 0.7, then it is certain that $\mathbf{s}$ is not $\theta$-similar to $\mathbf{q}$ because $0.6<0.7$. The verification algorithm can stop and avoid the rest of the scan. Note that when $\mathbf{q}$ is sparse, most of the time we have $lb(\tilde{\mathbf{s}})=\sum_{\tilde{\mathbf{s}}[i]\neq\bot}\tilde{\mathbf{s}}[i]\cdot\mathbf{q}[i]$ due to the existence of 0’s.

**Performance Guarantee. ** Next, we show that when the data is skewed, partial verification achieves a strong performance guarantee. We assume that each candidate $\mathbf{s}$ is stored in the database such that $\mathbf{s}[1]\geq\mathbf{s}[2]\geq\dots\geq\mathbf{s}[d]$ and we scan $\mathbf{s}$ sequentially starting from $\mathbf{s}[1]$. We notice that partial verification saves data accesses when there is a gap between the true similarity score and the threshold $\theta$. Intuitively, when this gap is close to 0, both the upper and lower bounds converge to $\theta$ as we scan $\mathbf{s}$ so we might not stop until the last element. When the gap is large and $\mathbf{s}$ is skewed, meaning that the first few values account for most of the candidate’s weight, then the first few $\mathbf{s}[i]\cdot\mathbf{q}[i]$ terms can provide large enough information for the lower/upper bounds to decide $\theta$-similarity. Formally,

Theorem B.4**.**

Suppose that a vector $\mathbf{s}$ is skewed: there exists an integer $k\leq d$ and constant value $c$ such that $\sum_{i=1}^{k}\mathbf{s}[i]^{2}\geq c$. For every query $(\mathbf{q},\theta)$, if $|\mathtt{cos}(\mathbf{s},\mathbf{q})-\theta|\geq\sqrt{1-c}$, then the number of accesses for verifying $\mathbf{s}$ is at most $k$.

Equivalently, Theorem B.4 says that if the true similarity is at least $\delta$ off the threshold $\theta$ (i.e. $|\mathtt{cos}(\mathbf{s},\mathbf{q})-\theta|\geq\delta$ for $\delta>0$), then it is only necessary to access $k$ entries of $\mathbf{s}$ with the smallest $k$ that satisfies $\sum_{i=1}^{k}\mathbf{s}[i]^{2}\geq 1-\delta^{2}$. For example, if $\delta=0.1$ and the first 20 entries of a candidate $\mathbf{s}$ account for >99% of $\sum_{i=1}^{d}\mathbf{s}[i]^{2}$, then it takes at most 20 accesses for verifying $\mathbf{s}$.

Proof B.5**.**

Case one: We first consider the case where $\mathtt{cos}(\mathbf{s},\mathbf{q})-\theta\geq\sqrt{1-c}$. In this case, we need to show $\sum_{i=1}^{k}\mathbf{s}[i]\cdot\mathbf{q}[i]\geq\theta$. Since $\mathtt{cos}(\mathbf{s},\mathbf{q})-\theta\geq\sqrt{1-c}$, which means

[TABLE]

it suffices to show that $\sum_{i=k+1}^{d}\mathbf{s}[i]\cdot\mathbf{q}[i]\leq\sqrt{1-c}$. This can be obtained by

[TABLE]

Case two:* We then consider the case where $\mathtt{cos}(\mathbf{s},\mathbf{q})-\theta\leq-\sqrt{1-c}$. In this case, we need to show *

[TABLE]

*The first term of the LHS is bounded by $\sum_{i=1}^{d}\mathbf{s}[i]\cdot\mathbf{q}[i]\leq\theta-\sqrt{1-c}$. The second term of the LHS is bounded by $\sqrt{1-\sum_{i=1}^{k}\mathbf{s}[i]^{2}}\leq\sqrt{1-c}$. Adding together, the LHS is bounded by $\theta$. *

Example B.6**.**

Figure 8a shows the number of accesses for each candidate during the verification phase of a real query with a vector having 100 non-zero valuesand $\theta=0.6$. The result shows that for most candidates, the number of accesses is much smaller than 100. In particular, 55.9% candidates need less than five accesses and 93.1% candidates need less than 30 accesses. This is because as shown in Figure 8b, only 0.23% of candidates have true similarity within $\pm$0.2 compared to $\theta$ (the range [0.4, 0.8]). The rest of the candidates, according to Theorem B.4, can be verified in a small number of steps.

Appendix C The baseline stopping condition is not tight

Theorem C.1**.**

The baseline stopping condition

[TABLE]

is complete but not tight.

Proof C.2**.**

For every position vector $\mathbf{b}$, $\varphi_{\mathsf{BL}}(\mathbf{b})=\mathtt{True}$ implies $\mathbf{q}\cdot L[\mathbf{b}]<\theta$. So for every $\mathbf{s}<L[\mathbf{b}]$, we also have $\mathbf{q}\cdot\mathbf{s}<\theta$ so $\varphi_{\mathsf{BL}}$ is complete.

To show the non-tightness, it is sufficient to show that for some position vector $\mathbf{b}$ where $\mathbf{b}$ is complete, $\varphi_{\mathsf{BL}}(\mathbf{b})$ is $\mathtt{False}$ so the traversal continues.

We illustrate a counterexample in Figure 9 with two dimensions (i.e., $d=2$). Given a query q, all possible $\theta$-similar vectors form a hyper-surface defining the set $\mathsf{ans}=\{\mathbf{s}|\ \|\mathbf{s}\|=1,\sum_{i=1}^{d}q_{i}\cdot s_{i}\geq\theta\}$. In Figure 9, $\|\mathbf{s}\|=1$ is the circular surface and $\sum_{i=1}^{d}q_{i}\cdot s_{i}\geq\theta$ is a half-plane so the set of points $\mathsf{ans}$ is the arc $\widehat{AB}$.

By definition, a position vector $\mathbf{b}$ is complete if the set $\{\mathbf{s}|\mathbf{s}<L[\mathbf{b}]\}$ contains no point in $\mathsf{ans}$. A position vector $\mathbf{b}$ satisfies $\varphi_{\mathsf{BL}}$ iff the point $L[\mathbf{b}]$ is above the plane $\sum_{i=1}^{d}q_{i}\cdot s_{i}=\theta$. It is clear from Figure 9 that if the point $L[\mathbf{b}]$ locates at the region BCD, then $\{\mathbf{s}|\mathbf{s}<L[\mathbf{b}]\}$ contains no point in $\widehat{AB}$ and is above the half-plane $\sum_{i=1}^{d}q_{i}\cdot s_{i}=\theta$. There exists a database of 2-d vectors such that $L[\mathbf{b}]$ resides in the BCD region for some position $\mathbf{b}$, so the stopping condition $\varphi_{\mathsf{BL}}$ is not tight.

Remark. The baseline stopping condition $\varphi_{\mathsf{BL}}$ is not tight because it does not take into account that all vectors in the database are unit vectors. In fact, one can show that $\varphi_{\mathsf{BL}}$ is tight and complete for inner product queries where the unit vector assumption is lifted. In addition, since $\varphi_{\mathsf{BL}}$ is not tight, any traversal strategy that works with $\varphi_{\mathsf{BL}}$ has no optimality guarantee in general since there can be a gap of arbitrary size between the stopping position by $\varphi_{\mathsf{BL}}$ and the one that is tight (i.e. there can be arbitrarily many points in the region $BCD$).

Appendix D Efficient computation of $\varphi_{\mathsf{TC}}$ with incremental maintenance

Testing the tight and complete condition $\varphi_{\mathsf{TC}}$ requires solving $\tau$ in Theorem 3.3, for which a direct application of the bisection method takes $\mathcal{O}(d)$ time. Next, we provide a more efficient algorithm based on incremental maintenance which takes only $\mathcal{O}(\log d)$ time for each test of $\varphi_{\mathsf{TC}}$.

According to the proof of Theorem 3.3,

[TABLE]

Wlog, suppose $\frac{L_{1}[b_{1}]}{q_{1}}\leq\dots\leq\frac{L_{d}[b_{d}]}{q_{d}}$ and $\tau$ is in the range $[\frac{L_{k}[b_{k}]}{q_{k}},\frac{L_{k+1}[b_{k+1}]}{q_{k+1}}]$ for some $k$. We have $s_{i}=L_{i}[b_{i}]$ for every $1\leq i\leq k$ and $s_{i}=q_{i}\cdot\tau$ for $k<i\leq d$. So if we let $\mathsf{eval}(k,\tau)$ be the function

[TABLE]

then for the largest $k$ such that $\mathsf{eval}(k,L_{k}[b_{k}]/q_{k})\leq 1$, $\tau$ can be computed by solving

[TABLE]

Then, $\mathsf{MS}(L[\mathbf{b}])$ can be computed as follows:

[TABLE]

Computing $\mathsf{MS}(L[\mathbf{b}])$ using the above approach directly requires that the $L_{i}[b_{i}]$’s are sorted in each step by $L_{i}[b_{i}]/q_{i}$, which requires $\mathcal{O}(d\log d)$ time. However, this is still too expensive as the stopping condition is checked in every step. Fortunately, we show that $\mathsf{MS}(L[\mathbf{b}])$ can be incrementally maintained in $\mathcal{O}(\log d)$ time as we describe below.

We use a binary search tree (BST) to maintain an order of the $L_{i}$’s sorted by $L_{i}[b_{i}]/q_{i}$. The BST supports the following two operations:

• •

$\mathsf{update(i)}$: update $L_{i}[b_{i}]\rightarrow L_{i}[b_{i}+1]$

• •

$\mathsf{compute()}$: return the value of $\mathsf{MS}(L[\mathbf{b}])$

The $\mathsf{compute()}$ operation essentially performs the binary search of finding the largest $k$ mentioned above. To ensure $\mathcal{O}(\log d)$ running time, we observe that from Equation (15) and (16), for any $k$, $\mathsf{MS}(L[\mathbf{b}])$ can be computed if $\sum_{i=1}^{k}L_{i}[b_{i}]\cdot q_{i}$, $\sum_{i=1}^{k}L_{i}[b_{i}]^{2}$, and $\sum_{i=1}^{k}q_{i}^{2}$ are available. Let $\mathsf{T}$ be the BST and each node in $\mathsf{T}$ is denoted as an integer $\mathsf{n}$, meaning that the node represents the list $L_{\mathsf{n}}$. We denote by $\mathsf{subtree(n)}$ the subtree of $\mathsf{T}$ rooted at node $\mathsf{n}$, and maintain the following information for each node $\mathsf{n}$:

• •

$\mathsf{n.key}$: $L_{n}[b_{n}]/q_{n}$,

• •

$\mathsf{n.LQ}$: $\sum_{i\in\mathsf{subtree(n)}}L_{i}[b_{i}]\cdot q_{i}$,

• •

$\mathsf{n.Q2}$: $\sum_{i\in\mathsf{subtree(n)}}q_{i}^{2}$ and

• •

$\mathsf{n.L2}$: $\sum_{i\in\mathsf{subtree(n)}}L_{i}[b_{i}]^{2}$ .

Thus, whenever there is a move on the list $L_{i}$, the key of the node $\mathsf{i}$ (i.e., $L_{i}[b_{i}]/q_{i}$) will be updated. Then we can remove the node $\mathsf{i}$ from the tree and insert it again using the new key, which takes $\mathcal{O}(\log d)$ time (with all the associated values being updated as well). To compute $\mathsf{MS}(L[\mathbf{b}])$, we need to handle two cases shown in Figure 10. In the first case, all the required information for computing $\mathsf{MS}(L[\mathbf{b}])$ are stored in the current node and its left subtree, while in the second case, the required information needs to be passed from the parent node to the current node to compute $\mathsf{MS}(L[\mathbf{b}])$.

Algorithm 2 shows the details to incrementally computes $\mathsf{MS}(L[\mathbf{b}])$ mentioned in Section D. For each node $\mathsf{n}$, we denote by $\mathsf{n.left}$ ($\mathsf{n.right}$) the left (right) child of $\mathsf{n}$.

Appendix E Traversal starting from the middle

Next, we consider a variant of the Gathering-Verification algorithm starting the traversal from the middle of the inverted lists instead of starting from the top like the classic $\mathsf{TA}$. In the most general setting, the starting position can be anywhere in the inverted lists, and two pointers (upper and lower) are used to indicate the traversed range for each inverted list. At each iteration, in addition to choosing the list to be traversed, the traversal strategy also decides whether the upper or the lower pointer of the selected list needs to be moved. We can show that in this setting, deciding a tight and complete stopping condition is np-hard thus intractable.

We start with some definitions. A configuration of the traversal is a pair of position vector $(\mathbf{a},\mathbf{b})$ where $\mathbf{a}$ and $\mathbf{b}$ are the upper and the lower position vector of the inverted lists $\{L_{i}\}_{i\in[d]}$ respectively. At each configuration, it is guaranteed that all vectors $\mathbf{s}$ satisfying $L_{i}[b_{i}]\leq s_{i}\leq L_{i}[a_{i}]$ for some dimension $i$ have been collected in the candidate set. So the tight stopping condition needs to test whether all the remaining vector $\mathbf{s}$ can be $\theta$-similar to the query $\mathbf{q}$. For each dimension $i$, the condition that $s_{i}$ needs to satisfy is

[TABLE]

which can be written into an equivalent quadratic form

[TABLE]

since $L_{i}[a_{i}]\geq L_{i}[b_{i}]$.

Thus, the tight stopping condition is equivalent to the following quadratic program:

[TABLE]

Let

[TABLE]

be the Lagrangian of (17) where $\lambda\in\mathbb{R}$ and $\mathbf{\mu}\in\mathbb{R}^{d}$ are the Lagrange multipliers. Then, the KKT conditions of (17) are

[TABLE]

in addition to the Primal feasibility in (17).

By the Complementary slackness condition, either $\mu_{i}=0$ or $(s_{i}-L_{i}[a_{i}])(s_{i}-L_{i}[b_{i}])=0$. If $\mu_{i}=0$, from the Stationarity condition, we have

[TABLE]

so $s_{i}=q_{i}/\lambda$.

If $\mu_{i}\neq 0$, then $(s_{i}-L_{i}[a_{i}])(s_{i}-L_{i}[b_{i}])=0$ so $s_{i}=L_{i}[a_{i}]$ or $s_{i}=L_{i}[b_{i}]$. Since $\mathbf{s}\cdot\mathbf{q}$ is maximized when $s_{i}$ is proportional to $q_{i}$ and $s_{i}\geq L_{i}[a_{i}]\lor s_{i}\leq L_{i}[b_{i}]$, the value of $s_{i}$ is determined as follows:

• •

if $q_{i}/\lambda\in[L_{i}[b_{i}],L_{i}[a_{i}]]$, then $s_{i}$ equals either the upper bound $L_{i}[a_{i}]$ or the lower bound $L_{i}[b_{i}]$;

• •

otherwise, $s_{i}=q_{i}/\lambda$.

Here, we can see that solving the quadratic program (17) is more difficult than the quadratic program (3) in Section 3 where the traversal starts from the top of the lists. Intuitively, when $s_{i}\neq q_{i}/\lambda$, unlike the solution for (3), there are still two possible options: the upper bound $L_{i}[a_{i}]$ and the lower bound $L_{i}[b_{i}]$, which lead to exponentially many combinations. One can easily show that even when $\lambda$ is fixed, computing the optimal $\mathbf{s}$ is np-hard by reduction from the subset sum problem [42]. This is because when $\lambda$ is fixed, the $s_{i}$’s that equal to $q_{i}/\lambda$ are fixed. When $L_{i}[b_{i}]=0$ for all $i$, checking whether $\sum_{i=1}^{d}s_{i}^{2}=1$ is essentially checking whether there exists a subset of $\{L_{i}[a_{i}]^{2}\}_{i\in[d]}$ whose sum is some fixed constant depending on $\lambda$.

Appendix F Proof of Theorem 4.5

Proof F.1**.**

Let $\{\mathbf{b}_{t}\}_{1\leq t\leq k}$ be the sequence of position vectors produced by the strategy $\mathcal{T}_{\mathsf{MR}}$.

Since each $\Delta_{i}$ is non-increasing and the strategy $\mathcal{T}_{\mathsf{MR}}$ chooses the dimension $i$ with the maximal $\Delta_{i}[b_{i}]$, then at each step $t$, the multiset $\{\Delta_{i}[j]|1\leq i\leq d,0\leq j\leq\mathbf{b}_{t}[i]\}$ contains the first $t$ largest values of all the $\Delta_{i}[j]$’s from the multiset $\{\Delta_{i}[j]|1\leq i\leq d,0\leq j<|L_{i}|\}$. Since the score $F(L[\mathbf{b}_{t}])$ equals

[TABLE]

it follows that for each $\mathbf{b}_{t}$ of $\mathcal{T}_{\mathsf{MR}}$, the score $F(L[\mathbf{b}_{t}])$ is the lowest score possible for any position vector reachable in $t$ steps. Thus, if the optimal access cost $\mathsf{OPT}$ is $t$ with an optimal stopping position $\mathbf{b}_{\mathsf{OPT}}$, then $\mathbf{b}_{t}$, the $t$-th position of $\mathcal{T}_{\mathsf{MR}}$, satisfies that $F(L[\mathbf{b}_{t}])\leq F(L[\mathbf{b}_{\mathsf{OPT}}])<\theta$. So $\mathcal{T}_{\mathsf{MR}}$ is optimal.

Appendix G Proof of Lemma 4.9

Proof G.1**.**

We construct a new collection of inverted lists $\{\tilde{L}_{i}\}_{1\leq i\leq d}$ from the original lists as follows. For every $i$ and every pair of consecutive indices $j_{k},j_{k+1}$ of $\mathsf{H}_{i}$, we assign

[TABLE]

Intuitively, we construct each $\tilde{L}_{i}$ by projecting the set of 2D points $\{(i,f_{i}(L_{i}[j]))\}_{j\in[j_{k},j_{k+1}]}$ onto the line passing through the two boundary points with index $j_{k}$ and $j_{k+1}$, which is essentially projecting the set of points onto the piecewise linear function defined by the convex hull vertices in $\mathsf{H}_{i}$ (See Figure 2(c) for an illustration). The new $\{\tilde{L}_{i}\}_{1\leq i\leq d}$ satisfies the following properties.

• (i)

By the construction of each convex hull $\mathsf{H}_{i}$, we have $\tilde{L}_{i}[j]\leq f_{i}(L_{i}[j])$ for every $i$ and $j$.

• (ii)

For every boundary position $\mathbf{b}$, we have $\tilde{L}[\mathbf{b}]=F(L[\mathbf{b}])$ since for every index $j$ on a convex hull $\mathsf{H}_{i}$, $\tilde{L}[j]=f_{i}(L_{i}[j])$.

• (iii)

The collection $\{\tilde{L}_{i}\}_{1\leq i\leq d}$ is ideally convex999where each $f_{i}$ of the decomposable function $F$ is the identity function*. In addition, the strategy $\mathcal{T}_{\mathsf{HL}}$ produces exactly the same sequence reduced by the max-reduction strategy $\mathcal{T}_{\mathsf{MR}}$ when $\{\tilde{L}_{i}\}_{1\leq i\leq d}$ is given as the input. By the same analysis for Theorem 4.5, for every position vector $\mathbf{b}$ generated by $\mathcal{T}_{\mathsf{HL}}$, $\tilde{L}[\mathbf{b}]$ is minimal among all position vectors reached within $|\mathbf{b}|$ steps.*

Combining (ii) and (iii), for every boundary position vector $\mathbf{b}$ generated by $\mathcal{T}_{\mathsf{HL}}$ and every $\mathbf{b}^{*}$ where $|\mathbf{b}^{*}|=\mathbf{b}$, we have $F(L[\mathbf{b}])=\tilde{L}[\mathbf{b}]\leq\tilde{L}[\mathbf{b}^{*}]$. Finally, by (i) and since $F$ is non-decreasing, we have $\tilde{L}[\mathbf{b}^{*}]\leq F(L[\mathbf{b}^{*}])$ so $F(L[\mathbf{b}])\leq F(L[\mathbf{b}^{*}])$ for every $\mathbf{b}^{*}$.

Appendix H Estimation of $\epsilon_{1}+\epsilon_{2}$

For $\epsilon_{1}$: A trivial upper bound of $\epsilon_{1}$ is $\tilde{\tau}-1$ since the initial value of $\tau$ is 1 and the gap $\tilde{F}(L[\mathbf{b}])-\mathsf{MS}(L[\mathbf{b}])$ is maximized when $\mathbf{b}=\mathbf{0}$. This upper bound can be improved as follows. We notice that in the proof of Theorem 4.12, given $|\mathbf{b}^{*}|=|\mathbf{b}_{l}|$, we need to have $F(L[\mathbf{b}^{*}])\geq\theta-\epsilon_{1}-\epsilon_{2}$ for every such $\mathbf{b}^{*}$, which is equivalent to requiring that this property holds for the $\mathbf{b}^{*}$ that minimizes $F(L[\mathbf{b}^{*}])$ given $|\mathbf{b}^{*}|=|\mathbf{b}_{l}|$. This $\mathbf{b}^{*}$ satisfies that $F(L[\mathbf{b}^{*}])\leq F(L[\mathbf{b}_{l}])$ and $F(L[\mathbf{b}_{l}])$ is known when the query is executed. This upper bound of $F(L[\mathbf{b}^{*}])$ implies a lower bound of $\tau$ at position $\mathbf{b}^{*}$, which also implies the following lower bound of $\mathsf{MS}(L[\mathbf{b}^{*}])-\tilde{F}(L[\mathbf{b}^{*}])$:

Lemma H.1**.**

Let $\mathbf{b}$ be an arbitrary position vector and let $\mathbf{b}^{*}$ be the position vector such that

[TABLE]

Then

[TABLE]

where $\tilde{F}$ is the decomposable function where each component is $f_{i}(x)=\min\{\tilde{\tau}\cdot q_{i},x\}\cdot q_{i}$ for constant $\tilde{\tau}$ and every $1\leq i\leq d$.

Proof H.2**.**

Let $\tau^{*}$ be the value of $\tau$ at $L[\mathbf{b}^{*}]$. We consider two cases separately: $\tau^{*}\geq\tilde{\tau}$ and $\tau^{*}<\tilde{\tau}$.

Case One: When $\tau^{*}\geq\tilde{\tau}$, since each $f_{i}$ increases as $\tilde{\tau}$ increases, we have

[TABLE]

thus $\mathsf{MS}(L[\mathbf{b}^{*}])\geq\tilde{F}(L[\mathbf{b}^{*}])$.

Case Two: Suppose $\tau^{*}<\tilde{\tau}$. We show the following Lemma.

Lemma H.3**.**

For every position $\mathbf{b}$ and constant $c$, $\mathsf{MS}(L[\mathbf{b}])<c$ implies that the $\tau$ at $L[\mathbf{b}]$ is at least $1/c$.

*Recall the notations $\mathsf{LQ}=\sum_{i:L_{i}[b_{i}]<\tau\cdot q_{i}}L_{i}[b_{i}]\cdot q_{i}$, $\mathsf{Q2}=\sum_{i:L_{i}[b_{i}]<\tau\cdot q_{i}}q_{i}^{2}$ and

$\mathsf{L2}=\sum_{i:L_{i}[b_{i}]<\tau\cdot q_{i}}L_{i}[b_{i}]^{2}$. Then $\tau$ satisfies that*

[TABLE]

and

[TABLE]

By rewriting Equation (19), we have $(1-\mathsf{Q2})=(1-\mathsf{L2})/\tau^{2}$. Plug this into (18), we have $\mathsf{LQ}+(1-\mathsf{L2})/\tau<c$. Since $L_{i}[b_{i}]<\tau\cdot q_{i}$, we have $\mathsf{LQ}>\mathsf{L2}/\tau$ so $1/\tau<c$ which means $\tau>1/c$. This completes the proof of Lemma H.3.

We know that since $\mathbf{b}^{*}$ minimizes $\mathsf{MS}(L[\mathbf{b}^{\prime}])$ among all $\mathbf{b}^{\prime}$ with $|\mathbf{b}^{\prime}|=|\mathbf{b}|$, we have $\mathsf{MS}(L[\mathbf{b}^{*}])\leq\mathsf{MS}(L[\mathbf{b}])$. By Lemma H.3, we have $\tau^{*}\geq 1/\mathsf{MS}(L[\mathbf{b}])$.

Now consider $\mathsf{MS}(L[\mathbf{b}^{*}])-\tilde{F}(L[\mathbf{b}^{*}])$. Since $\tau^{*}<\tilde{\tau}$, $\mathsf{MS}(L[\mathbf{b}^{*}])-\tilde{F}(L[\mathbf{b}^{*}])$ can be written as the sum of the following 3 terms:

• •

$\sum_{i:L_{i}[b_{i}]/q_{i}<\tau^{*}}(L_{i}[b_{i}]\cdot q_{i}-L_{i}[b_{i}]\cdot q_{i})$, which is always 0,

• •

$\sum_{i:\tau^{*}\leq L_{i}[b_{i}]/q_{i}<\tilde{\tau}}(q_{i}^{2}\cdot\tau^{*}-L_{i}[b_{i}]\cdot q_{i})$* and*

• •

$\sum_{i:\tilde{\tau}\leq L_{i}[b_{i}]/q_{i}}(q_{i}^{2}\cdot\tau^{*}-q_{i}^{2}\cdot\tilde{\tau})$.

In the second term, since each $L_{i}[b_{i}]$ is at most $q_{i}\cdot\tilde{\tau}$, so this term is greater than or equal to $\sum_{i:\tau^{*}\leq L_{i}[b_{i}]/q_{i}<\tilde{\tau}}(q_{i}^{2}\cdot\tau^{*}-q_{i}^{2}\cdot\tilde{\tau})$. Adding them together, we have $\mathsf{MS}(L[\mathbf{b}^{*}])-\tilde{F}(L[\mathbf{b}^{*}])$ is lower bounded by

[TABLE]

which is greater than or equal to $\tau^{*}-\tilde{\tau}$. Since $\tau^{*}\geq 1/\mathsf{MS}(L[\mathbf{b}])$, we have $\mathsf{MS}(L[\mathbf{b}^{*}])-\tilde{F}(L[\mathbf{b}^{*}])\geq 1/\mathsf{MS}(L[\mathbf{b}])-\tilde{\tau}$. Finally, combining case one and two together completes the proof of Lemma H.1.

Thus, when the query is given, $\epsilon_{1}$ is at most $\max\{0,\tilde{\tau}-1/\mathsf{MS}(L[\mathbf{b}_{l}])\}$.

For $\epsilon_{2}$: In general, there is no upper bound for $\tau$ since it can be as large as $L_{i}[b_{i}]/q_{i}$ for some $i$. The gap $\mathsf{MS}(L[\mathbf{b}])-\tilde{F}(L[\mathbf{b}])$ can be close to 1. However, the proof of Theorem 4.12 only requires $\epsilon_{2}$ to be at least the difference between $\mathsf{MS}$ and $\tilde{F}$ at position $\mathbf{b}_{l}$. So $\epsilon_{2}$ is upper-bounded by $\mathsf{MS}(L[\mathbf{b}_{l}])-\tilde{F}(L[\mathbf{b}_{l}])$ for a given query.

Summarizing the above analysis, the approximation factor $\epsilon$ is determined by the following two factors:

how much the approximation $\tilde{F}$ is smaller than the scoring function $\mathsf{MS}$ at $\mathbf{b}_{l}$, the starting position of the last segment chosen during the traversal, and 2. 2.

how much $\tilde{F}$ is bigger than $\mathsf{MS}$ at the optimal position with exactly $|\mathbf{b}_{l}|$ steps that minimizes $\mathsf{MS}$.

This yields the following upper bound of $\epsilon$:

[TABLE]

Appendix I Experiments on Document and Image datasets

In addition to the mass spectrometry dataset presented in Section 4, we also verify the near-convexity assumption in a document dataset101010https://www.kaggle.com/ycalisar/hotel-reviews-dataset-enriched. The documents are 515k hotel reviews from the website booking.com. We apply the doc2vec111111https://radimrehurek.com/gensim/models/doc2vec.html machine learning model to convert each hotel review to a vector representation and build the inverted index and convex hulls from the generated vectors. Similar to the mass spectrometry dataset, we observe that the values in the inverted lists have the similar “near-convex” shape as shown in Figure 11.

In document datasets, a cosine threshold query can be interpreted as “finding similar documents of a given one”. We ran 100 randomly chosen queries on a subset of 10,000 reviews with threshold $\theta=0.6$. The total number of accesses is 4,762,040 and the total sizes of the last gap is 374,521 which means that the cost additional to the optimal access is no more than 7.86% of the overall access cost.

Next, we repeat the same experiment on an image dataset121212http://vis-www.cs.umass.edu/lfw/ containing 13,000 images of human faces collected from the web. We use the img2vec model131313https://github.com/christiansafka/img2vec to convert each image to a vector. Once again, we observe the same near-convex shape patterns in the inverted lists constructed from the vectors as shown in Figure 12.

In image datasets, a cosine threshold query can be used for finding similar images to a given image. With the same setting above, we ran 100 random queries of faces in the whole dataset with the same threshold $\theta=0.6$. The total access cost is 118,795,452, the total size of last gaps is 473,999, so the additional access cost compared to $\mathsf{OPT}$ is no more than 0.40% of the overall access cost.

Appendix J Applying the proposed techniques to Top-$k$ cosine queries

We show that the proposed stopping condition and traversal strategy can also be applied to top-k cosine queries. In the top-k setting, each query is a pair $(\mathbf{q},k)$ where $\mathbf{q}$ is a query vector and $k$ is an integer parameter. A query $(\mathbf{q},k)$ asks for the top $k$ database vectors with the highest cosine similarity with the query vector $\mathbf{q}$.

The proposed techniques can be applied by adapting the classic structure of $\mathsf{TA}$. The algorithm traverses the 1-d inverted lists and keeps track of the $k$-th highest similarity score among the gathered vectors. Note that to guarantee tightness, the exact similarity scores need to be computed online, which results in additional computation cost since the exact computation can be avoided for threshold queries using partial verification (Section B). The stopping condition can be adapted as follows.

Theorem J.1**.**

At a position vector $\mathbf{b}$, let $\theta_{k}$ be the $k$-th highest score among vectors on or above $\mathbf{b}$. Then the following stopping condition is tight and complete:

[TABLE]

The score $\mathsf{MS}(L[\mathbf{b}])$ can be computed using the same $\mathcal{O}(\log d)$ incremental maintenance algorithm in Section D. The lower bound $\theta_{k}$ needs to be updated when a new candidate is gathered. Computing the similarity score takes $\mathcal{O}(d)$ time and updating $\theta_{k}$ can be done in $\mathcal{O}(\log k)$ using a binary heap.

The hull-based traversal strategy for inner product threshold queries can be directly applied to top-k inner product queries. The following near-optimality result can be shown.

Theorem J.2**.**

For a top-k inner product query $(\mathbf{q},k)$, the access cost of the hull-based traversal strategy $\mathcal{T}_{\mathsf{HL}}$ on a near-convex database $\mathcal{D}$ is at most $\mathsf{OPT}+c$ where $c$ is the convexity constant of $\mathcal{D}$.

The hull-based strategy can also be applied to top-k cosine queries. Recall that for cosine threshold queries, the hull-based traversal strategy operates on the convex hulls of a decomposable approximation $\tilde{F}$ of $\mathsf{MS}$ where $\tilde{F}(L[\mathbf{b}])=\sum_{i=1}^{d}\min\{q_{i}\cdot\tilde{\tau},L_{i}[b_{i}]\}\cdot q_{i}$. As discussed in Section 4.4, the choice of the constant $\tilde{\tau}$ is ideally the value $\tau$ of the optimal traversal path, which is dependent on the threshold $\theta$. Therefore, for top-k cosine queries where the final threshold $\theta_{k}$ is unknown in advance, a bad choice of $\tilde{\tau}$ can lead to poor approximation. In a practical implementation, the constant $\tilde{\tau}$ can be made query-dependent and more carefully tuned. We leave these aspects for future work.

Appendix K Additional Related Work

In this section, we provide additional related work.

K.1 Approximate approaches

Besides LSH, there are many other approximate approaches for high-dimensional similarity search, e.g., graph-based methods [32, 25, 75], product quantization [39, 7, 49], randomized KD-trees [64], priority search k-means tree[56], rank cover tree [34], randomized algorithms [73], HD-index [8], and clustering-based methods[50, 19].

K.2 TA-family algorithms

Bast et al. studied top-k ranked query processing with a different goal of minimizing the overall cost by scheduling the best sorted accesses and random accesses [9]. However, when the number of dimensions is high, it requires a large amount of online computations to frequently solve a Knapsack problem, which can be slow in practice. Note that, [9] only evaluated the number of sorted and random accesses in the experiments instead of the wall-clock time. Besides, it does not provide any theoretical guarantee, which is also applicable to  [40]. Akbarinia et al. proposed BPA to improve TA, but the optimality ratio is the same as TA in the worst case [3]. Deshpande et al. solved a special case of top-k problem by assuming that the attributes are drawn from a small value space [23]. Zhang et al. developed an algorithm targeting for a large number of lists [79] by merging lists into groups and then apply $\mathsf{TA}$. However, it requires the ranking function to be distributive that does not hold for the cosine similarity function. Yu et al. solved top-k query processing in subspaces [78] that are not applicable to general cosine threshold queries.

Some works considered non-monotonic ranking functions [80, 76]. For example, [80] focused on the combination of a boolean condition with a regular ranking function and [76] assumed the ranking functions are lower-bounded. The cosine function has not been considered in this line of work.

K.3 CSS in Hamming Space

[27] and [58] studied cosine similarity search in the Hamming space where each vector contains only binary values. They proposed to store the binary vectors efficiently into multiple hash tables. However, the techniques proposed there cannot be applied since transforming real-valued vectors to binary vectors loses information, and thus correctness is not guaranteed.

K.4 Inner product search

The cosine threshold querying is also related to inner product search where vectors may not be unit vectors, otherwise, inner product search is equivalent to cosine similarity search.

Teflioudi et al. proposed the LEMP framework [68, 67] to solve the inner product search where each vector may not be normalized. The main idea is to partition the vectors into buckets according to vector lengths and then apply an existing cosine similarity search (CSS) algorithm to each bucket. This work provides an efficient way for CSS that can be integrated to the LEMP framework.

Li et al. developed an algorithm FEXIPRO [51] for inner product search in recommender systems where vectors might contain negative values. Since we focus on non-negative values in this work, the proposed techniques (such as length-based filtering and monotonicity reduction) are not directly applicable.

There are also tree-based indexing and techniques for inner product search [60, 20]. But they are not scalable to high dimensions [43]. Another line of research is to leverage machine learning to solve the problem [57, 31, 63]. However, they do not provide accurate answers. Besides, they do not have any theoretical guarantee.

K.5 Dimensionality reduction

Since many techniques are affected negatively by the large number of dimensions, one potential solution is to apply dimensionality reduction (e.g. PCA, Johnson-Lindenstrauss) [54, 53, 41] before any search algorithm. However, this does not help much if the dimensions are not correlated and there are no hidden latent variables to be uncovered by dimensionality reduction. For example, in the mass spectrometry domain, each dimension represents a chemical compound/element and there is no physics justifying a correlation. We applied dimensionality reduction to the data vectors and turned out that only 4.3% of dimensions can be removed in order to preserve 99% of the distance.

K.6 Keyword search

This work is different from keyword search although the similarity function is (a variant of) the cosine function and the main technique is inverted index [55]. There are two main differences: (1) keyword search generally involves a few query terms [70]; (2) keyword search tends to return Web pages that contain all the query terms [55]. As a result, search engines (e.g., Apache Lucene) mainly sort the inverted lists by document IDs [14, 24, 16] to facilitate boolean intersection. Thus, those algorithms cannot be directly applied to cosine threshold queries. Although there are cases where the inverted lists are sorted by document frequency (or score in general) that are similar to this work, they usually follow TA [24]. This work significantly improves TA for cosine threshold queries.

K.7 Mass spectrometry search

For mass spectrometry search, the state-of-the-art approach is to partition the spectrum database into buckets based on the spectrum mass (i.e., molecular weight) and only search the buckets that have the similar mass to the query [44]. Each bucket is scanned linearly to obtain the results. The proposed techniques in this work can be applied to each bucket and improve the performance dramatically.

Library search, wherein a spectrum is compared to a series of reference spectra to find the most similar match, providing a putative identification for the query spectrum [77]. Methods in library search often use cosine [47] or cosine-derived [71] functions to compute spectrum similarity. Many of the state-of-the-art library search algorithms, e.g., SpectraST [47], Pepitome [21], X!Hunter [18], and [72] find candidates by doing a linear scan of peptides in the library that are within a given parent mass tolerance and computes a distance function against relevant spectra. This is likely because the libraries that are searched against have been small, but spectral libraries are getting larger, e.g., MassIVE-KB now has more than 2.1 million precursors. This work fundamentally improves on these, as it uses inverted lists with many optimizations to significantly reduce the candidate spectra from comparison. M-SPLIT used both a prefiltering and branch and bound technique to reduce the candidate set for each search [69]. This work improves on both these by computing the similar spectra directly, while still applying the filtering.

Database search [28] is where the experimental spectrum is compared to a host of theoretical spectra that are generated from a series of sequences. The theoretical spectra differ from experimental spectra in that they do not have intensities for each peak and contain all peaks which could be in the given peptide fragmentation. In other words, the values in the vectors are not important and the similarity function evaluates the number of shared dimensions that are irrelevant to the values. While the problem is different in practice than that described in this paper, it is still interesting to consider optimizations. Dutta and Chen proposed an LSH-based technique which embeds both the experimental spectra and the theoretical spectra onto a higher dimensional space [26]. While this method performs well at the expense of low recall rate, our method is guaranteed to not lose any potential matches. There are also other methods optimizing this problem by using preindexing [65, 44] but none do this kind of preindexing for calculating matches between experimental spectra.

K.8 Others

This work is also different from  [52] because that work focused on similarity search based on set-oriented p-norm similarity, which is different from cosine similarity.

In the literature, there are skyline-based approaches [48] to solve a special case of top-k ranked query processing where no any ranking function is specified. However, those approaches cannot be used to solve the problems that only involves unit vectors such that all the vectors belong to skyline points.