Computational Geometry with Probabilistically Noisy Primitive Operations

There is considerable prior work on sorting and searching with noisy comparison errors. For example, Feige, Raghavan, Peleg, and Upfal [15] show that one can sort in $O(n\log n)$ time w.h.p. with probabilistically noisy comparisons. Dereniowski, Lukasiewicz, and Uznanski [9] study noisy binary searching, deriving time bounds for constant factors involved. A similar study has also been done by Wang, Ghaddar, Zhu, and Wang [41]. Klein, Penninger, Sohler, and Woodruff solve linear programming in 2D under a different model of primitive errors [25], in which errors persist regardless of whether primitives are recomputed.

Other than the work by Groz et al. [19] mentioned above, we are not aware of prior work on computational geometry algorithms with random, independent noisy primitives. Nevertheless, considerable prior work has designed algorithms that can deal with input degeneracies, data imprecision, and arithmetic round-off errors. For example, several researchers have studied general methods for dealing with degeneracies in inputs to geometric algorithms, e.g., see [42, 13]. Researchers have designed algorithms for geometric objects with imprecise positions, e.g., see [27, 28]. In addition, significant prior work has dealt with arithmetic round-off errors and/or performing geometric primitive operations using exact arithmetic, e.g., see [16, 30]. While these prior works have made important contributions to algorithm implementation in computational geometry, they are orthogonal to the probabilistic noise model we consider in this paper.

Emamjomeh-Zadeh, Kempe, and Singhal [14] and Dereniowski, Tiegel, Uznański, and Wolleb-Graf [10] explore a generalization of noisy binary search to graphs, where one vertex in an undirected, positively weighted graph is a target. Their algorithm identifies the target by adaptively querying vertices. A query to a node $v$ either determines that $v$ is the target or produces an edge out of $v$ that lies on a shortest path from $v$ to the target. As in our model, the response to each query is wrong independently with probability . This problem is different than the graph search we study in this paper, however, which is better suited to applications in computational geometry. For example, in computational geometry applications, there is typically a search path, $P$, that needs to be traversed to a target vertex, but the search path $P$ need not be a shortest path. Furthermore, in such applications, if one queries using a node, $v$, that is not on $P$, it may not even be possible to identify a node adjacent to $v$ that is closer to the target vertex.

Viola [40] uses a technique similar to ours to handle errors in a communication protocol. In one problem studied by Viola, two participants with $n$-bit values seek to determine which of their two values is largest. This can be done by a noisy binary search for the highest-order differing bit position. Each search step performs a noisy equality test on two prefixes of the inputs, by exchanging single-bit hash values. The result is an $O(\log n)$ bound on the randomized communication complexity of the problem. Viola uses similar protocols for other problems including testing whether the sum of participant values is above some threshold. The noisy binary search protocol used by Viola directs the participants down a decision tree, with an efficient method to test whether the protocol has navigated down the wrong path in order to backtrack. One can think of our main technical lemma as a generalization of this work to apply to any DAG.

This work centers around our novel technique, path-guided pushdown random walks, described in Section 3. It extends noisy binary search in two ways: it can handle searches where the decision structure of comparisons is in general a DAG, not a binary tree, and it also correctly returns a non-answer in the case that the query value is not found. These two traits allow us to implement various geometric algorithms in the noisy comparison setting.

However, to apply path-guided pushdown random walks, one must design an oracle that, given a sequence of comparisons, can determine if one of them is incorrect in constant time. Because different geometric algorithms use different data structures and have different underlying geometry, we must develop a unique oracle for each one. The remainder of the paper describes noise-tolerant implementations with optimal running times for plane-sweep algorithms, point location, convex hulls in 2D and 3D, and Delaunay triangulations. We also present a dynamic 2D hull construction that runs in $O({\log }^{2}n)$ time w.h.p. per operation. See Table 1 for a summary of our results.

Our algorithms to construct the trapezoidal decomposition, Delaunay triangulation, and 3D convex hull are adaptions of classic randomized incremental constructions. A recent paper by Gudmundsson and Seybold [20] published at SODA 2022 claimed that such constructions produce search structures of depth $O(\log n)$ and size $O(n)$ w.h.p. in $n$, implying that the algorithms run in $O(n\log n)$ time w.h.p. Unfortunately, the authors have contacted us privately to say that their analysis of the size of such structures contains a bug. In this paper, we continue to use the standard expected time analysis.

We begin by noting that we can implement binary search trees storing geometric objects to support searches and updates in $O(\log n)$ time w.h.p. by adapting the noisy balanced binary search trees recently developed for numbers in quantum applications by Khadiev, Savelyev, Ziatdinov, and Melnikov [24] to our noisy-primitive framework for geometric objects. Plane-sweep algorithms require us to store an ordered sequence of events that are visited by the sweep-line throughout the course of the algorithm while maintaining an active set of geometric objects for the sweep line in a balanced binary tree. Some plane-sweep algorithms, such as the one of Lee and Preparata that divides a polygon into monotone pieces [26], have a static set of events that can simply be maintained as a sorted array that is constructed using an existing noisy sorting algorithm [15]. Others, however, add or change events to the event queue as the algorithm progresses. In addition, the algorithms must also maintain already-processed data in a way that allows for new data swept over to be efficiently incorporated. The two most common dynamic data structures in plane-sweep algorithms are dynamic binary search trees and priority queues. Throughout this section, noise is associated with the comparison “$a\le b$?” where $a$ and $b$ are geometric objects.

Once a data structure is built, its underlying structure can be manipulated without noise. For example, we need no knowledge of the values held in a tree to recognize that it is imbalanced and to perform a rotation, allowing implementation of self-balancing binary search trees. For example, Khadiev, Savelyev, Ziatdinov, and Melnikov [24] show how to implement red-black trees [21] in the noisy comparison model for numbers. We observe here that the same method can be used for ordered geometric objects compared with noisy geometric primitives, with a binary search tree serving as the search DAG and a root-to-leaf search path as the path. See Appendix B for an explanation for how to instantiate path-guided pushdown random walks on a binary search tree.

We can implement an event queue as a balanced binary search tree with a pointer to the smallest element to perform priority queue operations in the noisy setting.

Here, we present another plane-sweep application in the noisy setting.

We begin by observing that it is easy to construct two-dimensional convex hulls in $O(n\log n)$ time w.h.p. Namely, we simply sort the points using a noise-tolerant sorting algorithm [15] and then run Graham scan [8]. The Graham scan phase of the algorithm uses $O(n)$ calls to primitives, so we can simply use the trivial repetition strategy from Section 2 to implement this phase in $O(n\log n)$ time w.h.p. We note that this is the best possible, however, since just computing the minimum or maximum of $n$ values has an $\mathrm{\Omega }(n\log n)$-time lower bound for a high-probability bound in the noisy primitive model [15].

Overmars and van Leeuwen [32] show how to maintain a dynamic 2D convex hull with $O({\log }^{2}n)$ insertion and deletion times, $O(\log n)$ query time, and $O(n)$ space. They maintain a binary search tree, $T^{\ast }$, of points stored in the leaves ordered by $y$-coordinates and augment each internal node with a representation of half the convex hull of the points in that subtree. They define an $lc$-hull of $P$ as the convex hull of $P\cup \{(\mathrm{\infty },0)\}$, which is the left half of the hull of $P$ ($rc$-hull is defined symmetrically). It is easy to see that performing the same query on both the $lc$- and $rc$-hull of $P$ allows us to determine if a query point is inside, outside, or on the complete convex hull. We follow this same approach, showing how their approach can be implemented in the noisy-primitive model using path-guided pushdown random walks.

Overmars and van Leeuwen show that two $lc$-hulls $A$ and $C$ separated by a horizontal line can be merged in $O(\log n)$ time. Rather than being a single binary search, however, their method involves a joint binary search on the two trees that represent $A$ and $C$ with the aim of finding a tangent line that joins them. At each step, we follow one of ten cases to determine how to navigate the trees of $A$ and/or $C$. Because each decision may only advance in one tree, noisy binary search is not sufficient. Fortunately, as we show, the path-guided pushdown random walk framework is powerful enough to solve this problem in the noisy model. The challenge, of course, is to define an appropriate transition oracle. Say that the transition oracle has the ability to identify which case corresponds to our current position in $A$ and $C$. This decision determines which direction to navigate in $A$ and/or $C$ and corresponds to choosing a child to navigate down an implicit decision DAG wherein each node has ten children. It remains to show how a transition oracle acting truthfully can determine if we are on the correct path.

We have each node of $A$ and $C$ maintain $v.l$ and $v.r$ pointers as described in Appendix B. The values at $v.l$ and $v.r$ are ancestors of $v$ whose keys are just smaller and just larger than $v$ respectively. They determine the interval of possible values of any node in $v$’s subtree, which corresponds to an interval of points on each hull. See Figure 3 for an example.

Say we are currently at node $a\in A$ and node $c\in C$. We can perform four case comparisons: $a.r$ with $c.r$, $a.r$ with $c.l$, $a.l$ with $c.r$, and $a.l$ with $c.l$. Each outputs a region of $A$ and $C$ that the two tangent points must be located in. If the intersection of all four case comparisons contains the regions bounded by $a.r$ and $a.l$ as well as $c.r$ and $c.l$, then we are on a valid path by the correctness of Overmars and van Leeuwen’s case analysis [32]. Since the noise-free process takes $O(\log n)$ time, with a path-guided pushdown random walk using the transition oracle described above, this process takes $O(\log n)$ time w.h.p.

To query a hull, we can perform noisy binary search on the $lc$- and $rc$-hull structures using the transition oracle of Appendix B. Updating the convex hull utilizes the technical lemma discussed above along with split and join operations on binary search trees [32]. We conclude the following.

In this section, we show how to construct 3D convex hulls in $O(n\log n)$ time w.h.p. even with noisy primitive operations. The main challenges in this case are first to define an appropriate algorithm in the noise-free model and then define a good transition oracle for such an algorithm. For example, it does not seem possible to efficiently implement the divide-and-conquer algorithm of Preparata and Hong [35] in the noisy primitive model as each combine step performs $O(n)$ primitive operations. We begin by constructing a tetrahedron through any four points; we need the orientation of these four points, which can be found by the trivial repetition strategy in $O(\log n)$ time. The points within this tetrahedron can be discarded, again using the trivial repetition strategy in $O(\log n)$ time per point.

The main challenge to designing a transition oracle for the history DAG in this method is its nested nature. We have the nodes in the DAG themselves representing facets and pointing to the facets that replaced them when they were deleted. We also have the radial search structures, which maintain the set of new facets produced in a given iteration in sorted clockwise order around the point that generated them. See Figure 4 for a depiction of this.

Our transition oracle needs to be able to determine if we are on the correct path. Say we are given a query point, $q_i$, being inserted in the $i$th iteration of the RIC algorithm and are at node $v$ of the history DAG. To see if we are performing the right radial search, we first check if the ray from the origin to $q$ passes through $v$’s parent (in the history DAG, not the radial search structure). Then we must determine if we are performing radial search correctly. To do so, we can augment the radial search tree as described in Appendix B and perform two more orientation tests to determine if $v$ is on the right path in our current radial search structure. If a comparison returns false, then we use our stack to backtrack up the radial search structure or, if we are at its root, back to the parent node in the history DAG and its respective radial search tree. Thus, in just three comparisons, we can determine if we are on the right path in the history DAG. Similar to our trapezoidal decomposition analysis, we observe that the point-location cost in each iteration is $O({\mathrm{\ell }}_i+\log n)$, where ${\mathrm{\ell }}_i$ is the length of the unique path in the history DAG corresponding to our query $q$. An analysis by Mulmuley [31] shows that, over all $n$ iterations, ${\sum }_{i=1}^n{\mathrm{\ell }}_i=O(n\log n)$ in expectation. As a result, total point-location cost is $O(n\log n)$ in expectation.

Once a conflicting facet is found through the random walk, we can again walk around the hull and perform the trivial repetition strategy of Section 2 to determine the set of all conflicting facets, $X_i$, in $O(|X_i|\log n)$ time. Each facet must be created before it is destroyed, so ${\sum }_{i=1}^n|X_i|$ is upper-bounded by the total size of the history DAG. The same analysis by Mulmuley [31] shows that the expected size of the history DAG is $O(n)$. We conclude the following.

In this section, we show the following:

A shape Delaunay tessellation is a generalization of a Delaunay triangulation first defined by Drysdale [12]. Rather than a circle, we instantiate some convex compact set $C$ in ${𝐑}^2$ and redefine the empty circle property on an edge $\overline{pq}$ like so: given a set of points $S$ and some shape $C$, the edge $\overline{pq}$ exists in the shape Delaunay tessellation $D{T}_C(S)$ iff there exists some homothet of $C$ with $p$ and $q$ on its boundary that contains no other points of $S$ [4]. We will consider smooth shapes, as arbitrary non-smooth shapes may cause $D{T}_C(S)$ to not be a triangulation under certain point sets [4]. Note that all $L_p$ metrics for correspond to smooth shapes (rather than a unit circle, they are unit rounded squares) and furthermore, all smooth shapes produce triangulations of the convex hull of $S$ [4, 39]. Similar to above, we use the general position assumption that no four points lie on a homothet of $C$.

After inserting a point, $p_r$, in the incremental algorithm, we draw new edges that connect $p_r$ to the points that comprise the triangle(s) that $p$ landed in. These new edges are Delaunay as either they are fully enclosed by the homothets of $C$ that covered the triangle(s) $p$ was located in or they are enclosed by the union of previously-empty shapes that covered each edge of the triangles. After this, the algorithm finds adjacent edges, tests if they obey the empty shape property, and flips them if not. In Theorem 3 of a work by Aurenhammer and Paulini [4], the authors prove that for any convex shape $C$, local Delaunay flips lead to $D{T}_C(S)$ by showing that there exists some lexicographical ordering to these triangulations similar to how the Delaunay triangulation in $L_2$ maximizes the minimum angle over all triangles. Thus, the flipping algorithm used after an incremental insertion will terminate at a triangulation that is $D{T}_C(P_r)$, where $P_r$ is the set $\{p_1,\mathrm{\dots },p_r\}$.

Because the flipping process works as before, the point location structure can be easily adapted to this setting. As stated above, using a smooth shape $C$ guarantees that $D{T}_C(S)$ is an actual triangulation as opposed to a tree. Thus, in every iteration, there exists a triangle that the new point is enclosed by. The history DAG data structure is agnostic to the specific edge flips being made during the course of the algorithm and so works as expected. Our transition oracle also behaves the same.

In the Euclidean case, the backwards analysis of [8] uses the fact that a Delaunay triangulation of $r$ points has $O(r)$ edges, any new triangle created in iteration $r$ is incident to the newly inserted point $p_r$, and triangles are defined by at most three points. By planarity and by construction of the algorithm, all three hold true in this generalized case. Thus, we have the desired result.