我有一个由3D点组成的数组：
a = np.array([[2., 3., 8.], [10., 4., 3.], [58., 3., 4.], [34., 2., 43.]])
如何计算这些点的几何中心？
3个回答
35
我实现了Yehuda Vardi和Cun-Hui Zhang的几何中位数算法，该算法在他们的论文“多元L1-中位数及其相关数据深度”中有描述。所有内容都使用numpy向量化处理，因此速度应该非常快。我没有实现权重 - 只有无权重点。
import numpy as np
from scipy.spatial.distance import cdist, euclidean
def geometric_median(X, eps=1e-5):
y = np.mean(X, 0)
while True:
D = cdist(X, [y])
nonzeros = (D != 0)[:, 0]
Dinv = 1 / D[nonzeros]
Dinvs = np.sum(Dinv)
W = Dinv / Dinvs
T = np.sum(W * X[nonzeros], 0)
num_zeros = len(X) - np.sum(nonzeros)
if num_zeros == 0:
y1 = T
elif num_zeros == len(X):
return y
else:
R = (T - y) * Dinvs
r = np.linalg.norm(R)
rinv = 0 if r == 0 else num_zeros/r
y1 = max(0, 1-rinv)*T + min(1, rinv)*y
if euclidean(y, y1) < eps:
return y1
y = y1
除了默认的SO许可条款外，如果您喜欢，我还可以将上述代码发布在zlib许可下。
- orlp
8
我目前正在寻找替代使用 ArcGIS 的方法，这个方法非常接近于使用他们的 Median Center Geoprocessing 工具来计算中心点（700多个点精度在1米内）。您能否评论一下代码中发生了什么？我对 Python 和 Numpy 不是很熟练。我已经输入了一个包含pts xy坐标的numpy数组 X。如果您有时间解释一下函数中发生了什么，我将不胜感激。谢谢。 - Clubdebambos
1@Clubdebambos 我甚至不确定你想让我解释什么。你想知道为什么会有+-1m的差异吗？我无法访问他们的代码，所以我不知道他们在做什么。即使我知道，我的代码和他们的代码都有很多原因（数值不稳定性、浮点误差、多个可能的中位数候选项、算法或实现中的错误），而且我真的懒得去找出到底是哪一个原因。 - orlp
哦，这样啊。抱歉，我认为逐行注释可能不是很有用，你可以阅读我提供的论文。 - orlp
@Nabla，我已经有一段时间没有写代码了，但我认为应该与问题中的格式相同。 - orlp
我对上面的代码有同样的问题。我正在尝试在我的项目中实现这个geometric_median。我有一个包含WorkerID、Latitude和Longitude的numpy数组。我想要获取每个工人的中位数纬度和经度。X应该是由纬度和经度组成的元组吗？[y]又是什么呢？ - salvationishere
显示剩余3条评论
9
在这个gist或者从OpenAlea软件（CeCILL-C许可证）复制的下面函数中，使用了Weiszfeld的迭代算法来计算几何中位数。
import numpy as np
import math
import warnings
def geometric_median(X, numIter = 200):
"""
Compute the geometric median of a point sample.
The geometric median coordinates will be expressed in the Spatial Image reference system (not in real world metrics).
We use the Weiszfeld's algorithm (http://en.wikipedia.org/wiki/Geometric_median)
:Parameters:
- `X` (list|np.array) - voxels coordinate (3xN matrix)
- `numIter` (int) - limit the length of the search for global optimum
:Return:
- np.array((x,y,z)): geometric median of the coordinates;
"""
# -- Initialising 'median' to the centroid
y = np.mean(X,1)
# -- If the init point is in the set of points, we shift it:
while (y[0] in X[0]) and (y[1] in X[1]) and (y[2] in X[2]):
y+=0.1
convergence=False # boolean testing the convergence toward a global optimum
dist=[] # list recording the distance evolution
# -- Minimizing the sum of the squares of the distances between each points in 'X' and the median.
i=0
while ( (not convergence) and (i < numIter) ):
num_x, num_y, num_z = 0.0, 0.0, 0.0
denum = 0.0
m = X.shape[1]
d = 0
for j in range(0,m):
div = math.sqrt( (X[0,j]-y[0])**2 + (X[1,j]-y[1])**2 + (X[2,j]-y[2])**2 )
num_x += X[0,j] / div
num_y += X[1,j] / div
num_z += X[2,j] / div
denum += 1./div
d += div**2 # distance (to the median) to miminize
dist.append(d) # update of the distance evolution
if denum == 0.:
warnings.warn( "Couldn't compute a geometric median, please check your data!" )
return [0,0,0]
y = [num_x/denum, num_y/denum, num_z/denum] # update to the new value of the median
if i > 3:
convergence=(abs(dist[i]-dist[i-2])<0.1) # we test the convergence over three steps for stability
#~ print abs(dist[i]-dist[i-2]), convergence
i += 1
if i == numIter:
raise ValueError( "The Weiszfeld's algoritm did not converged after"+str(numIter)+"iterations !!!!!!!!!" )
# -- When convergence or iterations limit is reached we assume that we found the median.
return np.array(y)
或者，您可以使用C实现，该实现在此答案中提到，并使用例如ctypes将其与Python进行接口。
- rth
9
这看起来很慢，因为它是纯Python - 我正在寻找一个快速的numpy/scipy解决方案。 - orlp
1@orlp，快速的numpy/scipy解决方案需要代码可以向量化。乍一看，这似乎根本不可能。我想这里的问题是：你需要多快的速度？我猜写一个
cython版本的要点已经给你带来了非常好的速度。使用此答案中建议的c实现甚至可能更快。 - cel@Aprillion相比于使用numpy进行向量化，Python速度非常慢。 - orlp
看起来这个函数对于[[1, 1, 0], [1, -1, 0], [-1, -1, 0], [-1, 1, 0]]返回了错误的答案。 - shoosh
显示剩余4条评论
3
问题可以使用
对于 Nelder-Mead 方法，您将看到以下内容。要最小化的函数是从所有点到距离函数，即：
对于@orlp的方法，
对于Nelder mead，
scipy中的minimize模块轻松地进行近似处理。在此模块中，它提供了各种优化算法，从Nelder-Mead到Newton-CG。如果您不想麻烦高阶导数，则Nelder-mead算法特别有用，但代价是失去了一些精度。尽管如此，您只需要知道要最小化的函数即可让nelder-mead algorithm工作。现在，参考问题中的相同数组，如果我们使用@orlp的方法，我们将得到以下结果：
geometric_median(a)
# array([12.58942481, 3.51573852, 7.28710661])
对于 Nelder-Mead 方法，您将看到以下内容。要最小化的函数是从所有点到距离函数，即：
这里是代码：
from scipy.optimize import minimize
x = [point[0] for point in a]
y = [point[1] for point in a]
z = [point[2] for point in a]
x0 = np.array([sum(x)/len(x),sum(y)/len(y), sum(z)/len(z)])
def dist_func(x0):
return sum(((np.full(len(x),x0[0])-x)**2+(np.full(len(x),x0[1])-y)**2+(np.full(len(x),x0[2])-z)**2)**(1/2))
res = minimize(dist_func, x0, method='nelder-mead', options={'xtol': 1e-8, 'disp': True})
res.x
# array([12.58942487, 3.51573846, 7.28710679])
请注意，我使用“均值”作为算法的初始值。结果与@orlp的方法非常接近，但更加精确。正如我所提到的，您会有所牺牲，但仍然可以获得相当好的近似值。
Nelder Mead算法的性能
为此，我生成了一个带有10000个从以3.2为中心的正态分布中抽取的点的test_array。因此，几何中位数应该非常接近于[3.2, 3.2, 3.2]。
np.random.seed(3)
test_array = np.array([[np.random.normal(3.2,20),
np.random.normal(3.2,20),
np.random.normal(3.2,20)] for i in np.arange(10000)])
对于@orlp的方法，
%timeit geometric_median(test_array)
# 12.1 ms ± 270 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
# array([2.95151061, 3.14098477, 3.01468281])
对于Nelder mead，
%timeit res.x
# 565 ms ± 14.8 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
# array([2.95150898, 3.14098468, 3.01468276])
@orlp的方法非常快，而Nelder mead也不错。然而，Nelder mead方法是通用的，而@orlp的方法是特定于几何中位数的。你想选择的方法将取决于你的目的。如果你只想要一个近似值，我会轻松选择Nelder。如果你想要精确的结果，@orlp的方法既更快又更准确。
- shiba
网页内容由stack overflow 提供, 点击上面的可以查看英文原文，
原文链接
原文链接
np.median的实现方式，它似乎是基于partition实现的 - 这对于几何中位数来说不可能是正确的。 - orlp