python - 一种快速计算非空区域的方法

Question

我正在编写一些代码，选择 n 个通过原点的 5 维随机超平面。然后，它在单位球体上均匀地随机采样 no_points 点，并计算超平面创建的区域中有多少至少有一个点。使用以下 Python 代码可以相对简单地完成此操作。

import numpy as np

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

n = 100
d = 5
hpoints = points_on_sphere(n, d).T
for no_points in xrange(0, 10000000,100000):
    test_points = points_on_sphere(no_points,d).T 
    #The next two lines count how many of the test_points are in different regions created by the hyperplanes
    signs = np.sign(np.inner(test_points, hpoints))
    print no_points, len(set(map(tuple,signs)))

不幸的是，我计算不同区域有多少点的天真方法很慢。总的来说，该方法需要 O(no_points * n * d) 时间，并且在实践中，一旦 no_points 达到大约 1000000，它就会太慢并且太耗 RAM .特别是它在 no_points = 900,000 时达到 4GB RAM。

是否可以更有效地完成此操作，以便 no_points 可以相当快地达到 10,000,000(实际上，如果它可以达到 10 倍，那就太好了)并且使用不到 4GB 的 RAM ？

Answer 1

存储每个测试点如何相对于每个超平面进行分类是大量数据。我建议对点标签进行隐式基数排序，例如，

import numpy as np


d = 5
n = 100
N = 100000
is_boundary = np.zeros(N, dtype=bool)
tpoints = np.random.normal(size=(N, d))
tperm = np.arange(N)
for i in range(n):
    hpoint = np.random.normal(size=d)
    region = np.cumsum(is_boundary) * 2 + (np.inner(hpoint, tpoints) < 0.0)[tperm]
    region_order = np.argsort(region)
    is_boundary[1:] = np.diff(region[region_order])
    tperm = tperm[region_order]
print(np.sum(is_boundary))

此代码保留测试点的排列 (tperm)，以便同一区域中的所有点都是连续的。 boundary 指示每个点是否在排列顺序上与前一个不同的区域。对于每个连续的超平面，我们划分每个现有区域并有效地丢弃空区域以避免存储 2^100 个区域。

实际上，由于您有如此多的点而超平面却如此之少，因此不存储这些点更有意义。以下代码使用二进制编码将区域签名打包为两个 double 值。

import numpy as np


d = 5
hpoints = np.random.normal(size=(100, d))
bits = np.zeros((2, 100))
bits[0, :50] = 2.0 ** np.arange(50)
bits[1, 50:] = 2.0 ** np.arange(50)
N = 100000
uniques = set()
for i in xrange(0, N, 1000):
    tpoints = np.random.normal(size=(1000, d))
    signatures = np.inner(np.inner(tpoints, hpoints) < 0.0, bits)
    uniques.update(map(tuple, signatures))
print(len(uniques))