背景
对于给定的点列表,我想检查它们是否位于(可能是低维)单形的内部.
我想用Python做到这一点.
问题
编辑：从本质上讲,我想(反复)回答这个问题,你是否在于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
有没有办法提高准确性？
>我试着象征着同情,但不得不打断计算.花了太长时间.
>我不想再增加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)永远不会收敛到机器精度.事实上,您的代码正在做的是找到系统的唯一最小二乘解决方案的正确数值近似.
有很多方法可以解决系统无法提供解决方案的原因.一种方法是recall那个
A linear system
Ax=bis consistent if and only if the rank ofAis equal to the rank of[A|b], the matrixAaugmented withbas a column.
您可以按如下方式在数字上估计排名：
# 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.同样这样做也是可能的,但似乎很慢：
AA_b.rref()
