如何在高维球体表面相对均匀分布n个点

如何在一个高维球面尽可能均匀的分布$n$个点?

几天前导师让我思考这么一个问题,我的第一感觉是,寻找球内接正多面体,但是自己没找到数值化的推导,也没找到推广到高维方法;接着我又想:有没有类似’熵’的指标,类比最大熵用来衡量空间中向量分布的均匀程度?很遗憾这个也没找到 : (

恰好我知乎关注大牛王赟也思考过这个问题,并且找到了二维圆面、三维球面的解决方案、以及原理解释。这里我参考他文章中对二维圆面、三维球面解决方案,尝试给出高维球面的解决方案。

如何量化这个问题?可以想到的一种思路是:令点之间最短距离最大化(类似于支持向量机)。这样定义出来的问题叫做Tammes problem, 是密铺问题的一个特例。不幸的是,密铺问题往往没有很优雅的解。另一种定义方法是:把各个点看成同种电荷,让整个系统的电势能最小。这个问题可以通过模拟电荷的运动来解决,但计算复杂度非常高,而且只能得到数值解。

王赟在$StackOverflow$上面找到了一种可靠地近似解法:Evenly distributing n points on a sphere, 这种生成的点阵被称为斐波那契网格($Fibonacci lattice $或 $Fibonacci grid$)。至于为什么叫菲波那契网格,目前可以简单地用「黄金分割比也出现在菲波那契数列中」来解释,在下文中你还会见到菲波那契数的出现。这篇文章(Measurement of areas on a sphere using Fibonacci and latitude–longitude lattices)说明,使用菲波那契网格测量球面上不规则图形的面积,与用经纬网格(并加权)相比,误差可以减小 40%。

方法中运用黄金分割比的原理分析过程主要依赖连分式(Continued fraction),主要依赖于不同连分式的观察分析,这里具体可以看王赟文章

二维圆面

二维圆面中给出的公式如下(其中$\phi$为黄金分割比$\sqrt5 + 1\over 2$):

先来看看效果:

让我们来分析一下生成公式:

这两个公式从几何角度很好理解:对于第$n$个点,该点到原点的距离$\sqrt n$,然后该点所在角度为在上一个点(第$n-1$个点)的基础上旋转$2\pi\phi$,如此简单便完成近似均匀分布。

为什么这种思路会$work$呢?可以这样来看,对于第$k$个点,离原点距离为$\sqrt k$,它位于圆周上,则这个圆面积为$\pi k$,同理对于第$k-1$个点,园面积为$\pi (k-1)$,故第$k$个点所在圆环面积为$\pi$。由此可以看出,每一个点所占圆环面积均为$\pi$,这种情况保证了宏观上的均匀性:每$\pi$面积上就有一个点。但是这样还不够,不能保证环与环之间的点混乱均匀,对此就要用到黄金分割比保证不同环之间点选择的混乱均匀性,这里的原理具体可以看王赟文章

这里附上我的$python$实现代码:

1
2
3
4
5
6
7
def Fibonacci_2D(samples, radius):
n, r = samples, radius
indices = arange(0, samples, dtype=float) + 0.5
r_pts = sqrt(indices * r / n)
theta = pi *(5**0.5 - 1) * indices # theta为2pi_n_phi
x, y = r_pts*cos(theta), r_pts*sin(theta)
return [x, y]

三维球面

三维球面中点生成公式:

先来看看效果:

这三条公式对$\phi$值非常敏感。哪怕$\phi$稍微偏离黄金分割一点点,做出图的效果就不好。

下面解释一下这三条公式:

  1. $z_n = (2n-1)/N -1$(这里将球心归到原点),说明各个点竖坐标成等差序列,相当于把球面切成相同厚度的$N$层,并且在每一层厚度中点处的表面上取一个点。由于球冠表面积公式$S = 2 \pi Rh$,所以这样子切出来的各层有一个性质:侧面积都相等。这个性质保证了点阵在宏观上的均匀性:每$4\pi/N$面积上就有一个点。
  2. 第二、三条公式用到黄金分割比保证不同环之间点选择的混乱均匀性,这里的原理具体可以看王赟文章

下面附上我实现三维球面密铺的$python$代码:

1
2
3
4
5
6
7
8
9
def Fibonacci_3D(samples, radius, orign = [0,0,0]):
n, r = samples, radius
# + 0.5 表示在厚度的中点处取点
indices = arange(0, samples, dtype=float) + 0.5
theta = pi *(5**0.5 - 1) * indices
z = 2 * indices * r / n - r
x, y = sqrt(r ** 2 - (z ** 2))*cos(theta), sqrt(r ** 2 - (z ** 2))*sin(theta)
x, y, z = x+orign[0], y+orign[1], z +orign[2]
return [x, y, z]

高维球面

高维球面这部分是自己联想推导,不敢保证正确性。

$m$x$n$维矩阵来代表最终的$n$个点(行代表维度、列代表样本是生成函数生成矩阵方便的原因)

看二维圆面、三维球面可以看出,斐波那契网格主要思想就是通过等分某一维坐标将超球体等分为体积相等n个环,在环厚度中心位置取点,每次取点旋转角度通过黄金分割比来确定。比如说对于三维球体,选中第三维厚度中点后,剩下的只是一个二维圆周,只需要确定角度就可以了。那么对于四维呢,确定一维坐标后,剩下的三维组成球体,那么如何确定剩下的点?我的想法是对于高维球体某一维坐标确定后得到子超球体,继续用高维球体密铺的方法处理子超球体,只不过只取相对位的那一个点(比如当前要铺第五个点,那么在子超球体中,只取第五个点)。具体方法是假设需要密铺$m$维球体,前$m-2$维都用等分维度的方法确定坐标,但是要注意当前维度的超球体半径需要通过总半径减去前几维同一列(同一列即同一个样本)的平方和然后开根号($for$循环中的$temp_r$,比如$x^2+y^2+z^2+w^2=$,3确定$z=1,w=1$之后,$x、y$相应变为$x^2+y^2=1$)。对于矩阵的列代表不同的$n$个样本,每一列(每个样本)均对应不同的相对位置(用$pts_ratio$表示该列(样本)的分位数)。对于最后两维,用三维的方法确定即可,不过这里半径依然满足上述规则。

下面附上我实现高维球面密铺的$python$代码:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
def Fibonacci_ND(dimension, samples, radius ):
r = radius
res = [[] for i in range(dimension)]
# + 0.5 表示在厚度的中点处取点 pts_ratio记录不同点的相对位置。
pts_ratio = (arange(0, samples, dtype=float) + 0.5) / samples
res[0] = 2 * pts_ratio * r -r
# temp_dim 用于记录当前维度确定的坐标的平方和
temp_dim = (2 * pts_ratio * r -r) ** 2
for i in range(1, dimension-2):
# temp_r为确定前几维后剩下的超球体半径
temp_r = sqrt(r ** 2 - temp_dim)
new = 2 * temp_r * pts_ratio - temp_r
res[i] = new
temp_dim += new ** 2
temp_r = sqrt(r ** 2 - temp_dim)
# theta 为 2_pi_n_phi
theta = pi * ( 5 ** 0.5 - 1) * (arange(0, samples, dtype=float) + 0.5)
temp_x, temp_y = temp_r * cos(theta), temp_r * sin(theta)
res[dimension - 2], res[dimension - 1] = temp_x, temp_y
return res

验证

因为高维球面的密铺没法通过作图来目测验证,这里我用到自己写的函数(不知科学与否)来大概推测。

具体方案如下:因为大家的半径相等,所以点乘可以衡量两个点之间的角度差异,具体为差越小两个点越靠近。在这里我每次选取随机点,计算出这个随机点到其他点点乘值的排序$list$,选择最大点和第五大点添加到输出队列,最后$print$输出队列的最大值和最小值(两个最值可能来自不同随机点)看差别,如果差别很小暂可认为,算法较为均匀。

下面为我的$python$代码:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
## 嵌套list转置函数 ##
def mat_trans(list):
list2 = []
row, col = len(list), len(list[0])
for i in range(col):
list2.append([y[i] for y in list])
return list2
## 验证函数 ##
def angle_count(hh, n_rand, samples):
res = []
mat = mat_trans(hh)
# n_rand为随机数的个数,随机选择一个点,然后将这个点与其他点的点乘输出为
# 排序序列temp_list。因为半径都相等,点乘可以估量角度差值(越大差值越小)。
for i in range(n_rand):
temp_list = []
rand = random.randint(0, samples-1)
for j in range(samples):
if j != rand:
temp_dot = [mat[rand][i] * mat[j][i] for i in range(len(mat[0]))]
temp_list.append(sum(temp_dot))
temp_list.sort()
# 将每个随机点第一大跟第五大点乘值加入
res.append(temp_list[-1])
res.append(temp_list[-5])
return res

samp = Fibonacci_ND(4, 1000, 10)
kankan = angle_count(samp, 500, 1000)
# 输出最大和最小(可能来源不同点随机)
print(min(kankan), max(kankan))

简单试验,对于三维球面,样本量为1000,$min$和$max$误差在$0.7\%$;对于四维球面,样本量为1000,$min$和$max$误差在$1.4\%$。这里认为分布较为均匀。

参考资料:

  1. 10560 怎样在球面上「均匀」排列许多点?(上)
  2. 10561 怎样在球面上「均匀」排列许多点?(下)
  3. Evenly distributing n points on a sphere
  4. Measurement of areas on a sphere using Fibonacci and latitude–longitude lattices
0%