CODEWITHSUNDEEP
algorithmmath
Simd
Simd

Reputation: 21223

How can you find points in regions of space quickly

Consider 5 dimensional space which has been partitioned by at least two planes which go through the origin. If there are n planes then the number of distinct pyramidal regions this creates is exactly n^4/24-n^3/4+(23 n^2)/24-(3 n)/4+1 as long as they are in "general position", a term I make specific below.

A hyperplane through the origin can be defined by a single vector starting at the origin which is orthogonal to the plane. The hyperplanes are in general position if no three of the vectors defining the hyperplanes are coplanar.

How can I find one point per pyramidal region efficiently in 5 dimensional space? That is to return k points in 5d space, each of which is in a different pyramidal region. The value k will be part of the input.

My naive attempt so far just picks points uniformly at random on the surface of a hypersphere using (in Python)

def points_on_sphere(dim, N, norm=numpy.random.normal):
    """
    http://en.wikipedia.org/wiki/N-sphere#Generating_random_points
    """
    normal_deviates = norm(size=(N, dim))
    radius = numpy.sqrt((normal_deviates ** 2).sum(axis=0))
    points = normal_deviates / radius
    return points

And then chooses a subset so that each one is in a different region. I can verify that the points are really all in different pyramidal regions using the following code. 

def all_in_distinct_pyramidal_regions(testpoints, hyperplanes):
    signs = numpy.sign(numpy.inner(testpoints, hyperplanes))
    return (len(set(map(tuple,signs)))==len(signs))

Here is a picture of 4 planes in 3d with 14 points, one in each region (because 5d was hard to draw), one in each of the pyramidal regions created by the planes.

enter image description here

Is there an algorithm that works in 5 dimensional space to find k points, each in different pyramidal regions and takes anything like O(nk) time?

I have noticed that the pyramidal regions are convex. Convex spaces are closed under addition and in particular it seems that if you knew the hyperplanes that formed the border of a region you could potentially find a point inside the region by averaging points on the hyperplanes. However, I am not sure how to turn this observation into a full efficient algorithm.

Upvotes: 2

Views: 324

Answers (2)

David Eisenstat
David Eisenstat

Reputation: 65458

If you're willing to dig up a Python library for linear programming, then here's an approach that's easy to implement and reasonably efficient (time O(poly(n) k)). It's similar to the one that Gene proposed.

Let v1, ..., vm in R^5 \ 0 be normal vectors for some of the planes. To find a point in the intersection of their positive halfspaces, solve the linear program

maximize z
subject to
v1 . x - z >= 0
...
vm . x - z >= 0
x in [-1, 1]^5
z in R,

in variables x, z, where . denotes dot product. If z > 0 in the optimal solution, then the point x is in the intersection of the positive halfspaces.

To find one point in each region, enumerate recursively all vectors of length n containing entries ±1, with the following pruning. We interpret the entries of the vector as choosing the positive or negative halfspace of each normal vector. If the LP with a subset of constraints corresponding to the current prefix is infeasible, then prune the recursion. In Pythonish pseudocode:

def find_points(normals, prefix=()):
    "solve the LP with normals [sign * normals[i] for i, sign in enumerate(prefix)]"
    if "the LP is infeasible": pass
    elif len(normals) == len(prefix): yield x
    else:
        for point in find_points(normals, prefix + (1,)):  # or yield from
            yield point
        for point in find_points(normals, prefix + (-1,)):
            yield point

Thanks to the pruning, the number of internal nodes of the recursion tree is at most n times the number of regions, so this is not at all exponential.

The bottleneck of this approach is almost certainly solving the linear programs. Since we add constraints one at a time, if your LP solver supports it, specify the dual simplex algorithm and reuse the previous dual solution as the starting basis.

Upvotes: 1

Cybercartel
Cybercartel

Reputation: 12592

You can try a quadkey. Translate the points to a binary and concatenate it. Then sort the points. You can verify the upper bounds with the mostsignificant bits. You can find a source code from Microsoft Bing Maps:http://msdn.microsoft.com/en-us/library/bb259689.aspx.

Upvotes: 0

Related Questions