Python线性方程-高斯消去法
目标
在给定一组点的情况下，我试图找到满足所提供的所有点的线性方程的系数。
例如，如果我想要找到线性方程(ax + by +c= z)：
3x + 2y + 2 = z我至少需要三个三维点：
(2, 2, 12)
(3, 4, 19)
(4, 5, 24)给定足够的点和坐标(x，y，z)，我应该能够使用高斯消元法找到(a，b，c)。
然而，我认为在特殊情况下，我在解决矩阵问题时遇到了问题。您可以在这里查看我对python实现的第一次尝试：https://gist.github.com/anonymous/8188272
让我们来看几个例子。
数据集1
使用以下“手工制作”点(x，y，z)：
(2, 2, 12)
(3, 4, 19)
(4, 5, 24)对以下矩阵执行LU分解：
[[ 2. 2. 1. 12.]
[ 3. 4. 1. 19.]
[ 4. 5. 1. 24.]]反向求解U矩阵：
[[ 4. 5. 1. 24. ]
[ 0. -0.5 0.5 0. ]
[ 0. 0. 0.5 1. ]]返回结果(a，b，c)：
[3.0, 2.0, 2.0]对，是这样!一切看起来都很好。
数据集2
使用以下“手工制作”点(x，y，z)：
(3, 4, 19)
(4, 5, 24)
(5, 6, 29)对以下矩阵执行LU分解：
[[ 3. 4. 1. 19.]
[ 4. 5. 1. 24.]
[ 5. 6. 1. 29.]]反向求解U矩阵：
[[ 5.00000000e+00 6.00000000e+00 1.00000000e+00 2.90000000e+01]
[ 0.00000000e+00 4.00000000e-01 4.00000000e-01 1.60000000e+00]
[ 0.00000000e+00 0.00000000e+00 4.44089210e-16 0.00000000e+00]]返回结果(a，b，c)：
[1.0, 4.0, 0.0]虽然从技术上讲，这是一个解决方案，但不是我想要的！
数据集3
使用以下“手工制作”点(x，y，z)：
(5, 6, 29)
(6, 7, 34)
(7, 8, 39)对以下矩阵执行LU分解：
[[ 5. 6. 1. 29.]
[ 6. 7. 1. 34.]
[ 7. 8. 1. 39.]]反向求解U矩阵：
[[ 7.00000000e+00 8.00000000e+00 1.00000000e+00 3.90000000e+01]
[ 0.00000000e+00 2.85714286e-01 2.85714286e-01 1.14285714e+00]
[ 0.00000000e+00 0.00000000e+00 0.00000000e+00 3.55271368e-15]]实现崩溃...
Thoughts
在数据集2和3中，最后一行和倒数第二行是“特殊的”。倒数第二行的"b“和"c”具有相同的值(这在我的特殊示例中是真的！)。不幸的是，我缺乏数学知识来判断它是正面还是反面。
当最后一行全为零，并且它上面的行具有相等的值时，是否有什么特殊情况需要我处理？
提前感谢！
回答 3
Stack Overflow用户
发布于 2013-12-31 05:58:18
是的，这是一个需要以不同方式处理的特殊情况。在情况2和3中，您有一个rank deficient matrix。一般来说，它可能意味着有无限多的解决方案，或者没有解决方案。
您可以通过检查通过堆叠这些3向量生成的矩阵的determinant来确定这些情况是否会发生。
>>> import numpy as np
>>> from scipy.linalg import det
>>> data1 = np.array([(2, 2, 12), (3, 4, 19), (4, 5, 24)])
>>> data2 = np.array([(3, 4, 19), (4, 5, 24), (5, 6, 29)])
>>> data3 = np.array([(5, 6, 29), (6, 7, 34), (7, 8, 39)])
>>> det(data1)
-1.9999999999999982
>>> det(data2)
5.551115123125788e-17
>>> det(data3)
8.881784197001213e-16示例1是一个满秩矩阵，它从几何上告诉您这3个点是linearly independent。
示例2和3使矩阵的行列式为零，这告诉您这些点是线性相关的。
Stack Overflow用户
发布于 2013-12-31 05:33:41
看看numpy.linalg，scipy.linalg和scipy.optimize。
例如，使用数据集1
(2, 2, 12)
(3, 4, 19)
(4, 5, 24)对于3x + 2y + 2 = z，请使用numpy.linalg.solve
>>> a = np.array([[2, 2, 1],
[3, 4, 1],
[4, 5, 1]])
>>> b = np.array([12, 19, 24])
>>> np.linalg.solve(a, b)
array([ 3., 2., 2.])例如，使用数据集2
(3, 4, 19)
(4, 5, 24)
(5, 6, 29)我得到..。
>>> a = np.array([[3, 4, 1], [4, 5, 1], [5, 6, 1]])
>>> b = np.array([19, 24, 29])
>>> np.linalg.solve(a, b)
array([ 0.73333333, 4.26666667, -0.26666667])这就是你要找的答案吗？这是一个有效答案，但由于a的等级为2，因此[2., 3., 1.]也是一个有效答案。见下文..。
例如，使用数据集3
(5, 6, 29)
(6, 7, 34)
(7, 8, 39)重复过程...
>>> a = np.array([[5, 6, 1], [6, 7, 1], [7, 8, 1]])
>>> b = np.array([29, 34, 39])
>>> np.linalg.solve(a, b)
LinAlgError: Singular matrix这意味着系数的determinant是零，因此矩阵的inverse是无限的，或者一般而言，system of equations有无限多个解或没有解。
>>> np.linalg.det(a)
0.0
>>> 1. / np.linalg.det(a)
inf请记住，您是通过假设x = inv(a)b来求解系统ax = b的，因此inv(a)必须存在并且是有限的。如果a是奇异的，那么inv(a)是无穷大的。
>>> np.linalg.inv(a)
LinAlgError: Singular matrix因此，尝试使用least squares method来找到最佳解决方案如何？
>>> np.linalg.lstsq(a,b)
(array([ 2., 3., 1.]), # solution "x"
array([], dtype=float64), # residuals, empty if rank > a.shape[0] or < a.shape[1]
2, # rank
array([ 1.61842911e+01, 2.62145599e-01, 2.17200830e-16])) # singular values of "a"因此，它找到了最好的解决方案[2., 3., 1.]，幸运的是，它实际上是针对您的条件的解决方案！残差被作为空返回，因为正如@wim所说，a是秩不足的，EG:不等于方阵a__的维数，或全秩。
Stack Overflow用户
发布于 2013-12-31 05:54:59
您要寻找的是包含所有3个点的平面。对于数据集2和3，您的解决方案似乎表现得很奇怪的原因是，这三个点是共线的，因此不存在唯一的解决方案(即存在包含任何给定线的无限数量的平面)。这反映在你的LU分解中，因为你的矩阵是秩2，所以你的LU分解出现了一个‘0’行。
假设你坚持寻找包含所有3个点的平面，你需要确保你的矩阵实际上是3阶的。如果是，那么你就有了一个解决方案。如果它是秩2，则任何包含公共直线的平面都是有效解。
注意:如果你试图找到4个点的公共平面，那么你可能会发现不存在这样的解决方案。
https://stackoverflow.com/questions/20847853
复制相似问题