查找一个点是否位于一组点的凸包内，而不计算包本身

Question

测试点 P 是否位于由一组点 X 形成的凸包内的最简单方法是什么？

我想要一种在高维空间（例如，最多 40 维）中工作的算法，该算法不会显式计算凸包本身。有什么想法吗？

Answer 1

该问题可以通过找到线性规划的可行点来解决。如果您对完整细节感兴趣，而不是仅仅将 LP 插入现有求解器，我建议您阅读 Boyd 和 Vandenberghe 关于凸优化的优秀书籍。

设A = (X[1] X[2] ... X[n])，即第一列为v1，第二列为v2 等。

解决以下 LP 问题，

minimize (over x): 1
s.t.     Ax = P
         x^T * [1] = 1
         x[i] >= 0  \forall i

其中

1. x^T 是 x 的转置
2. [1] 是 all- 1 个向量。

问题有解当且仅当该点位于凸包内。

Answer 2

当且仅当从该点到其他点的所有向量的方向小于其周围圆/球/超球面的一半时，该点位于其他点的凸包之外。

下面是两点情况的草图，蓝色的点在凸包（绿色）内部，红色的点在凸包外部：

对于红色的，存在圆的二等分，使得向量从该点到凸包上的点只相交半个圆。
对于蓝点，不可能找到这样的二等分。

Answer 3

您不必计算凸包本身，因为它在多维空间中看起来相当麻烦。凸包有

一个众所周知的属性：任意向量（点）v 点 [v1, v2, .., vn] 的凸包内部可以表示为 sum(ki*vi)，其中 0 0 ;= ki <= 1 且 sum(ki) = 1。相应地，凸包之外的任何点都不会具有这样的表示。

在 m 维空间中，这将为我们提供具有 n 个未知数的 m 个线性方程组。

编辑
我不确定一般情况下这个新问题的复杂性，但对于 m = 2 来说，它似乎是线性的。也许，在这方面有更多经验的人会纠正我。

Answer 4

我在 16 维上也遇到了同样的问题。由于即使 qhull 也无法正常工作，因为必须生成太多面，所以我通过测试开发了自己的方法，是否可以在新点和参考数据之间找到分离超平面（我称之为“HyperHull”；）） 。

寻找分离超平面的问题可以转化为凸二次规划问题（参见：SVM） 。我使用 cvxopt 在 python 中完成了此操作，代码不到 170 行（包括 I/O）。即使存在问题，该算法在任何维度上都无需修改即可工作，即维度越高，船体上的点数越高（请参阅：关于多面体中随机点的凸包）。由于船体不是明确构造的，而是仅检查点是否在内部，因此与快速船体等相比，该算法在更高维度上具有很大的优势。

该算法可以“自然地”并行化，并且加速应该等于处理器的数量。

Answer 5

虽然最初的帖子是三年前的，但也许这个答案仍然会有帮助。 Gilbert-Johnson-Keerthi (GJK) 算法找到两个凸多胞体之间的最短距离，每个凸多胞体被定义为一组生成器的凸包——值得注意的是，凸包本身不需要计算。在一种特殊情况下，即所询问的情况，其中一个多面体只是一个点。为什么不尝试使用GJK算法来计算P与点X的凸包之间的距离？如果该距离为 0，则 P 在 X 内部（或至少在其边界上）。 Octave/Matlab 中的 GJK 实现（称为 ClosestPointInConvexPolytopeGJK.m）以及支持代码可在 http://www.99main.com/~centore/MunsellAndKubelkaMunkToolbox/MunsellAndKubelkaMunkToolbox.html 。 GJK 算法的简单描述可以在第 4 节中找到。 2 篇论文，位于 http://www.99main.com/~centore/ColourSciencePapers /GJKinConstrainedLeastSquares.pdf 。我已经对 31 维空间中的一些非常小的集合 X 使用了 GJK 算法，并且得到了很好的结果。 GJK 的性能与其他人推荐的线性规划方法相比如何还不确定（尽管任何比较都会很有趣）。 GJK 方法确实避免了计算凸包或用线性不等式表示包，这两者都可能非常耗时。希望这个答案有帮助。

Answer 6

您是否愿意接受通常应该有效但不能保证有效的启发式答案？如果你是，那么你可以尝试这个随机的想法。

令 f(x) 为到 P 的距离的立方乘以 X 中的事物数量，减去到 X 中所有点的距离的立方之和。从随机的某个位置开始，并使用爬山算法来最大化f(x) 对于距离 P 很远的球体中的 x。除了退化情况外，如果 P 不在凸包中，则应该很有可能找到 P 位于其一侧的超平面的法线，X 中的所有内容都在 的另一边。

Answer 7

使用 scipy.optimize.minimize 测试点是否位于外壳空间中的文章。

基于用户1071136的回答。

如果计算凸包，它确实会快得多，所以我为想要这样做的人添加了几行。我从 graham 扫描（仅限 2D）切换到 scipy qhull 算法。

scipy.optimize.minimize 文档：
https://docs.scipy.org/doc/scipy/reference/ Optimize.nonlin.html

import numpy as np
import scipy.optimize
import matplotlib.pyplot as plt
from scipy.spatial import ConvexHull


def hull_test(P, X, use_hull=True, verbose=True, hull_tolerance=1e-5, return_hull=True):
    if use_hull:
        hull = ConvexHull(X)
        X = X[hull.vertices]

    n_points = len(X)

    def F(x, X, P):
        return np.linalg.norm( np.dot( x.T, X ) - P )

    bnds = [[0, None]]*n_points # coefficients for each point must be > 0
    cons = ( {'type': 'eq', 'fun': lambda x: np.sum(x)-1} ) # Sum of coefficients must equal 1
    x0 = np.ones((n_points,1))/n_points # starting coefficients

    result = scipy.optimize.minimize(F, x0, args=(X, P), bounds=bnds, constraints=cons)

    if result.fun < hull_tolerance:
        hull_result = True
    else:
        hull_result = False

    if verbose:
        print( '# boundary points:', n_points)
        print( 'x.T * X - P:', F(result.x,X,P) )
        if hull_result: 
            print( 'Point P is in the hull space of X')
        else: 
            print( 'Point P is NOT in the hull space of X')

    if return_hull:
        return hull_result, X
    else:
        return hull_result

对一些示例数据进行测试：

n_dim = 3
n_points = 20
np.random.seed(0)

P = np.random.random(size=(1,n_dim))
X = np.random.random(size=(n_points,n_dim))

_, X_hull = hull_test(P, X, use_hull=True, hull_tolerance=1e-5, return_hull=True)

输出：

# boundary points: 14
x.T * X - P: 2.13984259782e-06
Point P is in the hull space of X

可视化它：

rows = max(1,n_dim-1)
cols = rows
plt.figure(figsize=(rows*3,cols*3))
for row in range(rows):
    for col in range(row, cols):
        col += 1
        plt.subplot(cols,rows,row*rows+col)
        plt.scatter(P[:,row],P[:,col],label='P',s=300)
        plt.scatter(X[:,row],X[:,col],label='X',alpha=0.5)
        plt.scatter(X_hull[:,row],X_hull[:,col],label='X_hull')
        plt.xlabel('x{}'.format(row))
        plt.ylabel('x{}'.format(col))
plt.tight_layout()