python - 超定义线性系统中的精度问题

Question

背景

对于给定的点列表，我想检查它们是否位于(可能是低维)单纯形的内部。我想在 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]

Answer 1

您在上面发布的系统没有精确的解决方案。所以 np.dot(AA,x) - u) 永远不会收敛到机器精度。事实上，您的代码所做的是为您的系统的唯一最小二乘解找到正确的数值近似值。

有多种方法可以查看系统为何无法求解。一种方法是 recall那个

A linear system Ax=b is consistent if and only if the rank of A is equal to the rank of[A|b], the matrix A augmented with b as 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。在 sympy 中这样做也是可能的，但似乎很慢:

AA_b.rref()