Python 线性方程 - 高斯消元法

Question

目标

给定一组点，我试图找到满足所提供的所有点的线性方程的系数。

例如，如果我想找到线性方程（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]]

实现崩溃...

想法

在数据集 2 和 3 中，最后一行和倒数第二行是“特殊的”。倒数第二行的“b”和“c”具有相同的值（在我的特殊示例中是这样！）。不幸的是，我缺乏从中得出正面或反面的数学知识。

当最后一行全为零并且它上面的行具有相等的值时，我需要处理一些特殊情况吗？

提前致谢！

Answer 1

退房numpy.linalg，scipy.linalg和scipy.optimize。

例如使用数据集1 h3>

(2, 2, 12)
(3, 4, 19)
(4, 5, 24)

for 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 h3>

(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 h3>

(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 987654326 @矩阵为无限，或一般术语，system of equations 987654327 @拥有无限数量的解决方案或没有解决方案。

>>> np.linalg.det(a)
0.0
>>> 1. / np.linalg.det(a)
inf

记住您正在通过假设x = inv(a)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.]，幸运的是你，实际上是您条件的解决方案！残差作为空返回，因为AS @WIM表示，a 987654351 @是排名缺陷，，例如：不等于方形矩阵a或完整等级 em>。

Answer 2

是的，这是一种特殊情况，您需要以不同的方式处理。在案例 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 制作行列式为零的矩阵，这告诉您这些点是线性相关的。

Answer 3

您正在寻找的是包含所有 3 个点的平面。对于数据集 2 和 3，您的解决方案似乎表现异常的原因是这三个点是共线的，因此不存在唯一的解决方案（即存在无限数量的包含任何给定线的平面）。这反映在您的 LU 分解中，出现“零”行，因为您的矩阵排名为 2。

假设你坚持找到包含所有 3 个点的平面，你需要确保你的矩阵实际上是 3 阶的。如果是，那么你有一个解决方案。如果是 rank 2，那么任何包含公共线的平面都是有效的解决方案。

注意：如果你试图找到 4 个点的公共平面，那么你可能会发现不存在这样的解。