超定义线性系统中的精度问题
Accuracy problems in overdefined linear system
背景
对于给定的点列表，我想检查它们是否位于（可能是低维）单纯形的内部。
我想在 Python.
中执行此操作
问题
编辑： 本质上，我想（反复）回答这个问题，u 是否在 A 的图像中（作为决定问题，所以只是 yes/no).
首先我做 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) 永远不会收敛到机器精度。事实上，您的代码正在做的是为您的系统的唯一最小二乘解找到正确的数值近似值。
有多种方法可以查看为什么系统无法求解。一种方法是 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()
背景
对于给定的点列表，我想检查它们是否位于（可能是低维）单纯形的内部。 我想在 Python.
中执行此操作问题
编辑： 本质上，我想（反复）回答这个问题，u 是否在 A 的图像中（作为决定问题，所以只是 yes/no).
首先我做 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) 永远不会收敛到机器精度。事实上，您的代码正在做的是为您的系统的唯一最小二乘解找到正确的数值近似值。
有多种方法可以查看为什么系统无法求解。一种方法是 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。在 sympy 中这样做也是可能的，但似乎很慢：
AA_b.rref()