Space-Efficient Approximate Spherical Range Counting in High Dimensions

We juxtapose the above problem definition with a different notion of approximate counting, adapted for example by Aronov and Har-Peled [10], where one is interested in approximating the number of points inside a query range. Given a query range $R_q$, the data structure shall return a value $\nu$ such that $(1-\epsilon )|R_q\cap P|\le \nu \le (1+\epsilon )|R_q\cap P|$. Aronov and Har-Peled [10] show that one can reduce such a counting query to a sequence of “emptiness” queries, determining if the query range is empty or not, of length polynomial in $\log |P|$ and ${\epsilon }^{-1}$. In contrast, our definition asks for an exact aggregate over a set whose geometry is allowed to vary within the annulus $B(q,1+\epsilon )\setminus B(q,1)$. This definition has been used, e.g., by Arya and Mount [11], although they study solutions in the low-dimensional regime.

In terms of known solutions, approximate spherical range counting in high dimensions is not clearly distinguished from other variants of approximate range searching. Perhaps the most promising approach for related problems in high dimensions is locality sensitive hashing (LSH) (see [26, 20, 6, 8]), which maps the points to buckets with the promise that near (resp. far) points are more (resp. less) likely to be mapped together. It is not clear how one could use an LSH-based approach to solve the counting variant without paying at least the cost of exploring all points in $B(q,1)$. This suggests that, within LSH-based frameworks, counting is unlikely to admit a fundamentally more efficient solution than reporting. In this work, we focus on designing solutions better suited to the counting variant.

Our main contribution is to design a data structure for the approximate range counting problem in high-dimensional Euclidean spaces. For any constant approximation parameter $\epsilon \in (0,1)$, our data structure uses near-linear space, polynomial preprocessing time, and query time that is sublinear for any sublinear number of points in the ambiguity zone, i.e., for any sublinear $t_q=|\left(B(q,1+\epsilon )\setminus B(q,1)\right)\cap P|$ the query time is sublinear.

Our approach is reminiscent of classic results for the exact range searching problem in low-dimensional Euclidean spaces. In particular, a central notion of our approach is that of partition trees [32, 14]. We adapt well-studied relations between spanning trees with low stabbing numbers (which also appears as “crossing” in previous works) and partition trees (see, e.g., [15]) to our approximate setting. A key difference in our setting is that we use a stronger notion of stabbing; we say that a point $q$ $\epsilon$-stabs a set $S$ if $\exists x\in S:\Vert x-q\Vert \le 1$ and $\exists y\in S:\Vert y-q\Vert \ge 1+\epsilon$ (classic stabbing would correspond to the case where $\epsilon =0$). An important step in our preprocessing algorithm is finding a “light” edge, i.e., a pair of input points that is $\epsilon$-stabbed by a polynomial in $d$ ratio of the relevant queries. Then, using the multiplicative weight update (MWU) method in a similar manner as [38, 15, 39], we are able to show that we can compute a spanning tree with low $\epsilon$-stabbing number, implying that we can build a partition tree such that for any query, we only need to visit a bounded number of nodes. To traverse the partition tree, we need to be able to decide for each internal node if the corresponding subset is $\epsilon$-stabbed by the query. We can only do this approximately, meaning that we may also explore nodes that are almost $\epsilon$-stabbed, but the overall cost of this exploration is proportional to $t_q$.

Our data structure is the first to achieve space usage that is near-linear in $n$ and polynomial in $d$, and query time that is sublinear in $n$ and polynomial in $d$ for any sublinear $t_q$. In Table 1, we juxtapose our result with previous known results that even implicitly solve the same problem. While LSH is originally proposed for the approximate near neighbor problem, where the data structure is only required to report one witness point within distance $1+\epsilon$, Ahle et al. [2] proposed an LSH-based solution for the problem of reporting all points. Their data structure is adaptive; while building time and space can be that of an LSH table prepared for worst-case inputs (parameter $\rho$ in Table 1), query time is that of the best LSH parameters for the query/instance (parameter $b$). Moreover, LSH tuned for near-linear space [8] can solve the reporting problem with only a logarithmic overhead in space and query time, since we can amplify the probability that every near neighbor will eventually be encountered. Since both LSH-based methods can solve the reporting variant, they can also solve the counting problem. However, it is not clear whether one can modify such an approach to avoid looking into all near points. Although our result does not guarantee sublinear query time in the worst case, it advances prior space-efficient work in high dimensions by making the complexity depend on the number of points in the annulus rather than on all points inside the entire ball. Trivially, the bucketing method, when combined with a Johnson-Lindenstrauss dimension reduction (as in [20]), can solve the counting problem; after reducing the dimension, one can store the relevant information in grid cells near the dataset. Such methods that discretize the query space, while providing efficient query times, unavoidably have huge space requirements. Finally, when the dimension is constant, partition trees [14] solve the problem exactly, while BBD-trees [11] leverage improved query times for approximate solutions (although still exponential in $d$).

In addition to our worst-case guarantees, we present an algorithm inspired by learning theory. Specifically, this algorithm is query-driven, i.e., it receives a number of query samples and aims to construct the data structure that yields the best expected query runtime for the distribution of the queries. This way, we design a much faster preprocessing algorithm (from $n^{\mathrm{𝗉𝗈𝗅𝗒}(1/\epsilon )}$, the complexity falls to $n^{O(1)}$), which gives us similar bounds for the average case.

Our approach is motivated by the viewpoint of analyzing algorithms through the lens of statistical learning (e.g., [19, 12]). While this has been influential in optimization and mechanism design, to the best of our knowledge, a similar treatment for data structures is absent. Although related ideas appear in learned index structures [31], LSH trees whose cutting rules are optimized with respect to the given dataset [7], learning-augmented binary trees with optimized search assuming access to advice from a frequency estimation oracle [30], these works do not directly cast the data structure construction itself as a PAC learning problem. We also remark that the related notions of self-improving [3] and instance-optimal algorithms [1] address different guarantees. A simplified version of our result follows:

We believe this is a neat way to use tools from other areas to design efficient algorithms for difficult problems or problems with only existence guarantees, and it would be interesting to see what other data-structure problems might benefit from this viewpoint.

In Section 2 we present notation and useful known results and prove auxiliary lemmas that will be used later. In Section 3 we analyze a data structure that supports approximate stabbing queries. Subsequently, in Section 4 we prove that a spanning tree with sublinear $\epsilon$-stabbing number exists and can be computed in polynomial time, which we then use to prove that a partition tree with sublinear query time can be efficiently constructed in Section 5. In Section 6 we present our data-driven algorithm for building a partition tree with near-optimal expected visiting number. Missing proofs can be found in the full version.

We build $L=c\log n$ independent data structures as in Lemma 9. For any query point $q\in {ℝ}^d$, we query the data structures, and we report as follows. If we find a point $x$ such that $\Vert x-q\Vert \le 1+\epsilon$ and a point $y$ such that $\Vert y-q\Vert \ge 1$ we return “stabbed”. Otherwise, if all returned points are within distance $1+\epsilon$, we return “covered”. Else, we return “disjoint”.

If there is a pair $x,y$ in $P$ such that $\Vert x-q\Vert \le 1$ and $\Vert y-q\Vert \ge 1+\epsilon$ then the probability that all $L$ data structures fail to return a point $x^{\prime }$ such that $\Vert x^{\prime }-q\Vert \le 1+\epsilon$ or a point $y^{\prime }$ such that $\Vert y^{\prime }-q\Vert \ge 1$ is by union bound at most $2\cdot {0.05}^L\le n^{-c}$. Moreover, if all points of $P$ are within distance $1$ from $q$ (resp. at a distance of at least $1+\epsilon )$, then all returned points are within distance $1+\epsilon$ from $q$ (resp. at a distance of at least $1$), implying that the query algorithm correctly identifies “covered” (resp. “disjoint”). $◀$

If $d\ge n^{\rho }$, then the statement trivially holds, since $|Q|=m$. Thus, we assume . For any two points $x,y$, let $S_{x,y}$ be the set of points in ${ℝ}^d$ that $\epsilon$-stab $x,y$. We define a range space $ℛ=({ℝ}^d,\{S_{x,y}\mid x,y\in {ℝ}^d,x\ne y\})$. Notice that for any point $p\in {ℝ}^d$, determining if $p\in S_{x,y}$ reduces to evaluating the signs of the following four polynomials of degree $2$ in $d$ variables: ${\Vert p-x\Vert }^2-1$, ${\Vert p-x\Vert }^2-{(1+\epsilon )}^2$, ${\Vert p-y\Vert }^2-1$, ${\Vert p-y\Vert }^2-{(1+\epsilon )}^2$. This implies (see, e.g., [9][Theorem 8.3]) that the VC dimension of $ℛ$ is $O(d)$. Hence, by Theorem 3, for any and constant $c>0$, a random sample $N_{\delta }$ of $O\left(\frac{d}{\delta }\left(\log \left(\frac{1}{\delta }\right)+\log n\right)\right)$ points of $Q$ is a $\delta$-net for $ℛ$ with probability $1-n^{-c-1}$. Moreover, we can find such a $\delta$-net by randomly sampling from $Q$, using the assumed multiset sampling data structure of Lemma 5, in $\tilde{O}(\frac{d}{\delta }\log n)\cdot O(\log |U|)$ time.

By Theorem 7, there exists a terminal embedding $f$ to dimension $d^{\prime }=O({\epsilon }^{-2}\log |N_{\delta }|)$ such that $\forall p\in N_{\delta },q\in {ℝ}^d$, $\Vert f(p)-f(q)\Vert \in (1\pm \epsilon /20)\Vert p-q\Vert$. The time needed to compute $f(P)$ is in $O(n)\cdot \mathrm{poly}(|N_{\delta }|)$ for the probability of success to be at least $1-n^{-c-1}$.

Now let $x,y\in P$ be $\epsilon$-stabbed by some $q\in Q$, and further assume wlog that $\Vert x-q\Vert \le 1$ and $\Vert y-q\Vert \ge 1+\epsilon$. We obtain $\Vert f(x)-f(q)\Vert \le 1+\epsilon /20$ and $\Vert f(y)-f(q)\Vert \ge (1+\epsilon )(1-\epsilon /20)\ge 1+3\epsilon /4$. Hence, $f(x),f(y)$ are $(\epsilon /2,1+\epsilon /20)$-stabbed by $f(q)$. Moreover, by the triangle inequality, $\Vert f(x)-f(y)\Vert \ge \epsilon /2$. Now consider the grid of side-length $\epsilon /(4\sqrt{d^{\prime }})$ in ${ℝ}^{d^{\prime }}$ and let $U_B$ be the set of all grid cells that intersect with at least one ball of radius $1+\epsilon$ centered at a point of $f(N_{\delta })$. There are

such cells, where $|N_{\delta }|=O\left(\frac{d}{\delta }\left(\log \left(\frac{1}{\delta }\right)+\log n\right)\right)$. Now let $\delta =\frac{d}{n^{\rho }}$ where $\rho =\mathrm{\Theta }({\epsilon }^2/\log (1/\epsilon ))$ is chosen sufficiently small (in terms of the hidden constants) so that . By the pigeonhole principle, there is either a pair of points of $P$ mapped to the same cell of $U_B$, under $f$, or there is a pair of points of $P$ that are not mapped to $U_B$. In either case, the two points are not $(\epsilon /2,1+\epsilon /20)$-stabbed by any point of $f(N_{\delta })$. To see this, notice that if the two points are in the same cell, then they are within distance $\epsilon /4$ from each other, contradicting the fact that they are $(\epsilon /2,1+\epsilon /20)$-stabbed (which would mean they are separated by a distance of at least $\epsilon /2$). If they are not mapped in $U_B$ then neither of them is within distance $1+\epsilon /20$ from any point of $f(N_{\delta })$, also implying that they cannot be $(\epsilon /2,1+\epsilon /20)$-stabbed by any point of $f(N_{\delta })$. Since the two points are not $(\epsilon /2,1+\epsilon /20)$-stabbed by any point of $f(N_{\delta })$ in the mapped space, their pre-images (under $f$) are not $\epsilon$-stabbed by any point of $N_{\delta }$. Since $N_{\delta }$ is a $\delta$-net, it must be that at most $\delta m$ points of $Q$ $\epsilon$-stab those two points.

We now focus on the running time needed to find a pair of points as above. By Theorem 7, we can map points under $f$ in $\mathrm{poly}(n)$ time. To compute the two points, we essentially use the bucketing method (see, e.g., [20]) for balls centered at $f(N_{\delta })$. We first store all cells of $U_B$ in a dictionary data structure. Then, for each point $x\in f(P)$, we query the data structure to see if its cell $C_x$ is stored, i.e., $C_x\in U_B$. If $C_x\in U_B$ then we store $x$ in a linked list in the bucket of $C_x$. We can stop and output any two points stored in the same linked list, or any two points whose cells do not belong to $U_B$. Iterating over all points of $f(P)$ costs $\tilde{O}(dn)$ assuming that the dictionary is implemented using balanced binary trees. The overall probability of success follows from a union bound over failure events of the random sampling for computing $N_{\delta }$ and the failure event of computing the terminal embedding $f$. $◀$

Knowing that a good pair exists, and that we can find it, we repeat the process $n/2$ times, in order to build a good forest, in the sense that it has a small $\epsilon$-stabbing number. We use the MWU technique introduced in the context of trees with low stabbing number in [15].

We conclude with the main result of this section, the tractability of a spanning tree with a sublinear $\epsilon$-stabbing number. This follows from $O(\log n)$ iterative applications of Lemma 12.

We first apply a JL transformation [28] to reduce the dimension to $O({\epsilon }^{-2}\log n)$. By Theorem 6, and by applying a union bound over $n$ vectors, for any $c>0$ and any query $q\in {ℝ}^d$, with probability $1-1/n^c$, all distances from $P$ are approximately preserved within factors $1\pm \epsilon /10$. This initial step costs $O(dn)$ time. For the rest of the proof, we assume that all points lie in ${ℝ}^{O({\epsilon }^{-2}\log n)}$. Then, we assume that queries are snapped to the grid ${𝒢}_{\gamma }^d$, $\gamma =\epsilon /(10\sqrt{d})$. This induces an additive $\pm \epsilon /10$ error to the distances between $P$ and the query. Hence, points within distance $1$ from the query may be moved to a distance up to $1+\epsilon /5$, and points at a distance of at least $1+\epsilon$ may be moved to a distance of at least $1+4\epsilon /5$. Notice that $1+\epsilon /2\le \frac{1+4\epsilon /5}{1+\epsilon /5}$, hence it suffices to use a more refined error of $\epsilon /2$ for the rest of the construction and assume that the space is rescaled by a factor $1/(1+\epsilon /5)$.

We build a binary partition tree $𝒯$ on the input with error parameter $\epsilon /2$ using Corollary 16. The construction of $𝒯$ can be done in $n^{O({\epsilon }^{-2}\log (1/\epsilon ))}$ time and be successful for a fixed $q$, in that $\zeta (q,𝒯)\le n^{1-\mathrm{\Theta }({\epsilon }^2/\log (1/\epsilon ))}+t_q$, with probability $1-1/n^c$ for any constant $c>0$. For each internal node $v$, let $P_v$ be the set of points in the leaves of the subtree rooted at $v$. In each internal node $v$, we build a data structure ${𝒟}_v$ as in Theorem 10 on $P_v$ with error parameter $\epsilon /2$, and we store the cumulative weight $s_v$ of the points in $P_v$, i.e., $s_v={\sum }_{p\in P_v}{w}_p$. By Theorem 10, each construction of a data structure ${𝒟}_v$, can be done in $\tilde{O}(d|P_v|)$ time and be successful with probability $1-1/n^c$ for any constant $c>0$. By a union bound over all inner nodes, we claim that all data structures ${𝒟}_v$ are also successful for $q$ with probability $1-\mathrm{\Theta }\left(\frac{n\log n}{n^c}\right)$. By Theorem 10, the time and space needed to build data structures for the inner nodes is near-linear per level of the tree (since each level defines a partition).

To answer a query $q\in {ℝ}^d$, we traverse the tree and collect all cumulative weights $s_v$ in each internal node $v$, whose associated data structure ${𝒟}_v$ returns “covered”. We continue traversing the children of a node $v$ if and only if ${𝒟}_v$ returns “stabbed”. By Corollary 16, assuming that the tree has been constructed correctly for $q\in {ℝ}^d$, the number of nodes visited during the query algorithm for input $q\in {ℝ}^d$ is $\tilde{O}(n^{\alpha }+t_q)$, where $\alpha =1-\mathrm{\Theta }({\epsilon }^2/\log (1/\epsilon ))$. By Theorem 10, the time needed to query a data structure ${𝒟}_v$ in an internal node $v$ is $\tilde{O}({|P_v|}^{\beta })$, where $\beta =1-\mathrm{\Theta }({\epsilon }^2)$. At each level $i$ of $𝒯$, the number of nodes whose data structures are queried is upper bounded by both $2^i$ and $\zeta (q,𝒯)$. Hence, the running time at level $i$ for query $q$ is $\tilde{O}(T_i^q(n))$, where $T_i^q(n)$ is defined as follows:

because if $i\le \log (n^{\alpha }+t_q)$ then $\min \{{\mathrm{\hspace{0.17em}2}}^i,n^{\alpha }+t_q\}=2^i$ and if $i>\log (n^{\alpha }+t_q)$ then $\min \{{\mathrm{\hspace{0.17em}2}}^i,n^{\alpha }+t_q\}=n^{\alpha }+t_q$. Now, since ,

where $\varrho =1-\beta$. We observe that the query algorithm will not explore any descendants of an internal node $v$ if ${𝒟}_v$ returns “disjoint” when ${\min }_{x\in P_v}\Vert x-q\Vert >1$, or if ${𝒟}_v$ returns “covered” when ${\max }_{x\in P_v}\Vert x-q\Vert \le 1+\epsilon /2$, in which case it takes into account $s_v$. Thus, the answer is calculated by summing over cumulative weights $s_{v_1},\mathrm{\dots },s_{v_{\mathrm{\ell }}}$ of nodes corresponding to subsets $P_{v_1},\mathrm{\dots },P_{v_{\mathrm{\ell }}}$ such that $B(q,1)\cap P\subseteq {\bigcup }_{i\in [\mathrm{\ell }]}{P}_{v_i}\subseteq B(q,1+\epsilon /2)\cap P$. $◀$

By Lemma 18, with probability at least $1-\delta$, for any spanning tree $T$ of $P$, $m^{-1}\cdot {\sum }_{q\in S}\sigma (q,T)$ is in

Now, let $E(T)$ denote the edges of spanning tree $T$ and $\tilde{T}={\arg \min }_T\{m^{-1}\cdot {\sum }_{q\in S}\sigma (q,T)\}$. Notice that ${\sum }_{q\in S}\sigma (q,T)={\sum }_{q\in S}{\sum }_{e\in E(T)}\sigma (q,\{e\})={\sum }_{e\in E(T)}{\sum }_{q\in S}\sigma (q,\{e\})$ for any spanning tree $T$. Hence, we can compute $\tilde{T}$, by first constructing a complete weighted graph on $P$ where we assign weight ${\sum }_{q\in S}\sigma (q,\{e\})$ to each edge $e\in \left(\genfrac{}{}{0pt}{}{P}{2}\right)$, and then computing the minimum spanning tree on that graph, using Kruskal’s algorithm. The running time of this step is $O(n^{2}m)$, which corresponds to the step of building the graph.