>Background
对于给定的点列表，我想检查它们是否位于(可能是较低维)单纯形的内部。 我想在Python中做到这一点。
问题
编辑：从本质上讲，我想(反复)回答这个问题，u是否在于A的形象(作为决策问题，所以只是是/否)。
首先我做QR分解，然后求解系统，最后检查解决方案是否正确。
import scipy.linalg
import numpy as np
Q,R = np.linalg.qr(AA)
for u in points:
x = scipy.linalg.solve_triangular(R, np.dot(Q.T, u))
print(all(x <= 1-1e-6) and all(x >= 1e-6) and all(abs(np.dot(AA,x) - u) < 1e-6))
但是，我遇到了数值问题，准确性似乎太差了。我有一个点，它位于内部(根据先前的线性程序计算)，但上面的代码无法识别这一点。
矩阵的条件约为 100，形状为(36,35)，所以看起来并不那么可怕，但误差略高于1e-6
有没有办法提高准确性？
- 我用
sympy尝试了它，但不得不中断计算。花了太长时间。 - 我不想再提高门槛
1e-6。 - 由于系统被过度定义(即单纯形具有较低的维度)，我需要最后一次检查，我的解决方案是否正确。
数据
AA = np.array([
[1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 2, 2, 2, 2, 4, 4, 6, 6, 6, 6, 6, 6, 6, 8],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 2, 2, 2, 4, 6, 0, 0, 0, 4, 4, 0, 0, 0, 0, 0, 0, 0, 0, 2, 6],
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 4, 0, 0, 0, 4, 4, 4, 0, 0, 4, 0, 2, 0, 0, 0, 0, 0, 2, 2, 6, 0, 0],
[0, 0, 0, 0, 0, 2, 2, 2, 4, 4, 6, 6, 0, 0, 0, 0, 2, 0, 0, 0, 0, 2, 0, 0, 0, 0, 2, 0, 0, 6, 0, 2, 4, 0, 4],
[0, 0, 0, 4, 10, 0, 0, 6, 0, 0, 2, 2, 2, 0, 2, 2, 4, 2, 2, 0, 0, 2, 4, 2, 0, 0, 0, 0, 4, 6, 2, 0, 2, 0, 0],
[0, 2, 2, 0, 4, 4, 4, 2, 2, 4, 2, 4, 2, 2, 2, 2, 2, 0, 0, 0, 0, 6, 2, 2, 0, 4, 0, 4, 0, 0, 0, 2, 0, 2, 0],
[0, 0, 2, 0, 4, 0, 2, 6, 0, 0, 0, 0, 6, 2, 0, 0, 0, 4, 6, 6, 6, 0, 4, 4, 0, 2, 6, 8, 0, 0, 0, 4, 0, 0, 0],
[0, 0, 0, 0, 6, 8, 0, 0, 2, 0, 2, 4, 2, 6, 2, 4, 2, 0, 0, 2, 2, 0, 0, 2, 4, 0, 0, 2, 4, 2, 10, 2, 0, 12, 0],
[0, 0, 0, 0, 0, 2, 0, 6, 2, 2, 0, 0, 2, 4, 0, 0, 0, 2, 0, 2, 2, 4, 0, 0, 6, 0, 0, 0, 0, 0, 0, 0, 2, 0, 0],
[0, 2, 0, 0, 2, 0, 0, 0, 0, 0, 4, 4, 0, 0, 0, 4, 2, 4, 8, 0, 0, 2, 2, 0, 0, 6, 0, 4, 0, 0, 2, 0, 0, 0, 0],
[0, 0, 2, 6, 4, 0, 0, 0, 0, 0, 8, 0, 6, 0, 4, 10, 0, 2, 0, 0, 4, 2, 6, 2, 2, 0, 2, 0, 0, 8, 2, 0, 0, 2, 0],
[0, 0, 2, 0, 4, 8, 2, 14, 0, 0, 6, 0, 0, 0, 4, 8, 0, 4, 2, 4, 0, 0, 0, 0, 0, 2, 8, 0, 0, 0, 0, 0, 0, 8, 2],
[0, 0, 6, 6, 0, 0, 0, 0, 4, 0, 0, 0, 0, 4, 4, 0, 0, 2, 0, 0, 2, 2, 2, 0, 0, 0, 2, 0, 0, 8, 0, 4, 0, 2, 0],
[0, 6, 0, 2, 0, 2, 0, 2, 0, 0, 0, 10, 0, 0, 0, 0, 4, 0, 0, 2, 0, 0, 0, 2, 0, 2, 0, 2, 2, 2, 2, 4, 0, 2, 0],
[0, 0, 0, 6, 0, 0, 8, 0, 0, 0, 8, 0, 8, 0, 0, 0, 2, 4, 0, 0, 2, 0, 0, 0, 2, 0, 4, 0, 0, 2, 0, 10, 8, 0, 2],
[0, 2, 2, 0, 0, 0, 2, 0, 2, 2, 0, 0, 0, 0, 0, 2, 0, 0, 2, 0, 2, 0, 2, 6, 2, 0, 0, 2, 6, 2, 4, 0, 4, 0, 2],
[0, 10, 0, 0, 0, 2, 2, 4, 0, 2, 0, 0, 0, 8, 0, 6, 6, 4, 0, 4, 2, 2, 0, 2, 0, 4, 0, 0, 0, 2, 0, 0, 4, 0, 10],
[0, 4, 0, 0, 4, 2, 0, 0, 0, 0, 2, 2, 0, 0, 0, 0, 0, 0, 4, 4, 0, 8, 2, 12, 4, 8, 0, 0, 4, 0, 6, 2, 10, 0, 2],
[0, 0, 4, 8, 6, 2, 0, 0, 12, 0, 2, 2, 0, 0, 0, 0, 2, 4, 0, 2, 0, 2, 0, 2, 4, 0, 2, 2, 10, 0, 0, 0, 0, 2, 2],
[0, 2, 2, 2, 2, 0, 0, 0, 2, 0, 0, 0, 6, 6, 2, 0, 2, 2, 4, 2, 4, 4, 4, 2, 10, 6, 2, 0, 2, 0, 0, 0, 2, 2, 8],
[0, 4, 0, 4, 0, 0, 0, 0, 2, 16, 0, 2, 0, 0, 6, 0, 4, 0, 0, 4, 0, 2, 0, 0, 0, 0, 0, 4, 0, 0, 0, 0, 0, 2, 0],
[0, 0, 0, 0, 0, 0, 2, 2, 0, 0, 2, 0, 6, 0, 0, 0, 2, 0, 0, 0, 2, 0, 4, 0, 2, 4, 2, 2, 2, 0, 0, 0, 0, 2, 0],
[0, 0, 4, 6, 6, 2, 8, 0, 4, 0, 0, 2, 0, 2, 0, 2, 4, 2, 0, 0, 2, 2, 0, 6, 2, 12, 4, 0, 0, 0, 4, 0, 0, 4, 2],
[0, 0, 4, 0, 0, 0, 0, 2, 4, 2, 2, 0, 0, 0, 2, 2, 4, 4, 4, 4, 2, 4, 4, 0, 2, 0, 0, 6, 2, 2, 2, 6, 0, 0, 2],
[0, 0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 10, 2, 0, 8, 2, 4, 0, 0, 0, 0, 0, 0, 0, 4, 0, 0, 4, 2, 2, 0, 0, 2, 0, 2],
[0, 2, 0, 0, 0, 8, 4, 2, 4, 6, 0, 0, 0, 2, 0, 0, 2, 0, 8, 0, 0, 0, 0, 0, 2, 0, 2, 0, 12, 4, 0, 0, 2, 0, 4],
[0, 0, 0, 0, 0, 0, 8, 2, 8, 2, 2, 2, 2, 0, 0, 2, 2, 10, 0, 2, 0, 0, 2, 2, 0, 2, 4, 0, 0, 6, 4, 4, 2, 4, 0],
[0, 4, 2, 0, 0, 0, 0, 0, 0, 0, 0, 2, 2, 4, 6, 0, 0, 0, 0, 4, 0, 2, 6, 2, 0, 0, 0, 0, 2, 0, 4, 0, 4, 0, 0],
[0, 0, 2, 0, 0, 2, 0, 0, 0, 6, 0, 0, 0, 4, 2, 0, 0, 0, 0, 0, 8, 0, 4, 0, 0, 0, 2, 0, 0, 0, 2, 2, 0, 2, 2],
[0, 8, 4, 4, 0, 0, 12, 0, 6, 0, 6, 2, 0, 6, 4, 2, 2, 0, 4, 0, 0, 2, 2, 0, 4, 0, 0, 0, 2, 0, 8, 0, 0, 0, 0],
[0, 0, 8, 0, 0, 0, 0, 0, 0, 4, 0, 2, 6, 2, 0, 0, 0, 0, 0, 2, 0, 2, 0, 0, 2, 2, 0, 2, 0, 2, 0, 4, 0, 0, 0],
[0, 6, 10, 6, 6, 4, 0, 0, 0, 4, 0, 2, 0, 4, 4, 0, 0, 4, 0, 0, 10, 0, 0, 2, 0, 2, 0, 10, 0, 0, 0, 0, 0, 2, 0],
[0, 2, 0, 4, 0, 8, 2, 0, 2, 4, 0, 2, 2, 0, 0, 4, 2, 0, 2, 4, 2, 4, 4, 0, 0, 0, 2, 0, 0, 0, 0, 0, 0, 2, 0],
[0, 4, 0, 0, 0, 0, 2, 10, 0, 2, 4, 0, 2, 0, 6, 4, 2, 0, 4, 0, 6, 2, 0, 2, 0, 0, 10, 0, 0, 0, 0, 6, 0, 2, 0],
[0, 0, 2, 2, 0, 4, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 2, 0, 2, 0, 0, 2, 0, 2, 0, 0, 2, 2, 0, 0, 0, 0, 2, 0, 2]])
v = np.array([1, 1, 1, 1, 2, 2, 2, 2, 2, 1, 2, 2, 2, 1, 2, 2, 1, 2, 2, 2, 2, 2, 1, 2, 2, 2, 2, 2, 1, 1, 2, 1, 2, 2, 2, 1])
points = [v]
您在上面发布的系统没有确切的解决方案。因此，np.dot(AA,x) - u)永远不会收敛到机器精度。事实上，您的代码正在做的是找到系统唯一最小二乘解的正确数值近似。
有多种方法可以了解系统无法提供解决方案的原因。 一种方法是回忆一下
线性系统
Ax=b是一致的，当且仅当A秩等于[A|b]秩，矩阵A用b作为列进行扩充。
您可以按如下数字估计排名：
# reshape the RHS to a column so we can combine it with AA
b = points[0].reshape((36,1))
# append the column to form the augmented matrix
AA_b = np.hstack((AA, b))
print("AA rank: %d" % np.linalg.matrix_rank(AA))
print("[AA|b] rank: %d" % np.linalg.matrix_rank(AA_b))
这表明AA的等级是 35，但[A|b]的等级是 36，因此系统无法解决。
您还可以通过将增强矩阵放入缩减的行梯队形式来说服自己这是真的。 我能够在wxMaxima中非常快速地完成此操作，并验证增强矩阵是否等同于标识的RREF。 在 sympy 中这样做也是可能的，但似乎很慢：
AA_b.rref()