如何快速判断n维点是否处于n维单纯形内？

我拥有大规模数据集，分析流程中需判断每个点是否处于n维单纯形内，希望在不使用并行化的前提下实现高效计算。数据集维度不固定，因此方案需适配任意维度。为便于演示，以下采用2D示例（理论上数学逻辑适用于任意维度）：

重心坐标法 #

最初尝试使用重心坐标法，但实现结果不可靠：

import numpy as np
import matplotlib.pyplot as plt

def is_in_simplex(point, T_inv, simplex):
    first_n = np.matmul(
        T_inv, (point - simplex[-1])
    )
    last = 1 - np.sum(first_n)
    bary = np.concatenate((first_n, [last]))
    return np.all((bary <= 1) & (bary >= 0))

# setup

simplex = np.array([[0, 0,], [8, 8,], [10, 3]])
rng = np.random.default_rng()
test_points = rng.random((10, 2))*10

# Maths starts here

T = np.array(simplex[:-1] - simplex[-1]).T
T_inv = np.linalg.inv(T)
within = np.vectorize(is_in_simplex, excluded={1, 2})(test_points, T_inv, simplex)

# drawing

polygon = np.concatenate([simplex, np.array([simplex[0]])])
print()
plt.plot(*polygon.T)
plt.scatter(*test_points.T)
for i, p in enumerate(test_points, 0):
    print(f"{i}\t{p}\t{test_points[i]}\t{within[i]}")
    plt.annotate(i, p)

问题与耗时：输出结果存在错误，如点3、5、6本应处于单纯形内，但重心坐标计算错误；is_in_simplex()返回数组而非单个布尔值。该方法处理10、100、1000、10000点的耗时分别为0.0383、0.0487、0.0994、0.523秒。

线性规划法 #

尝试使用线性规划法，结果正确但速度远低于预期：

import numpy as np
from scipy.optimize import linprog
import time

def vectorizable(point, simplexT, coeffs):
    b = np.r_[point, np.ones(1)]
    lp = linprog(coeffs, A_eq = simplexT, b_eq = b)
    return lp.success

dims = 2
rng = np.random.default_rng()
test_points = rng.random((10, dims))*10
simplex = np.array([[0, 0,], [8, 8,], [10, 3]])
coeffs = np.zeros(len(simplex))
simplex_T = np.r_[simplex.T,np.ones((1,len(simplex)))]

start_time = time.time()
in_simplex = np.vectorize(vectorizable,
                        excluded={1, 2},
                        signature="(n) -> ()")(test_points, simplex_T, coeffs)

print(f"----- {time.time() - start_time} seconds -----")

polygon = np.concatenate([simplex, np.array([simplex[0]])])
print()
plt.plot(*polygon.T)
plt.scatter(*test_points.T)
for i, p in enumerate(test_points, 0):
    print(f"{i}\t{p}\t{in_simplex[i]}")
    plt.annotate(i, p)

问题与耗时：该方法处理10000点耗时4-8秒，维度提升后耗时增至数十秒甚至数分钟。

希望优化重心坐标法的实现，或加速线性规划法，尽量避免并行化。

内容的提问来源于stack exchange，提问作者P. H. Allus