Python中是否有高斯消元的标准解决方案？

Question

在 scipy/numpy/... 宇宙中是否有标准的矩阵高斯消元方法？

通过谷歌可以找到许多代码片段，但如果可能的话，我更倾向于使用“可信赖”的模块。

- flonk

---

你是在寻找特别的高斯消元方法，还是其他解线性方程组、求逆矩阵等方式？ - Fred Foo

дёҚпјҢжҲ‘еҸӘйңҖиҰҒй«ҳж–Ҝж¶Ҳе…ғгҖӮеҺҹеӣ жҳҜпјҢжҲ‘жңүдёҖдёӘз”ұNдёӘж–№зЁӢз»„жҲҗзҡ„зі»з»ҹпјҢ其秩дёәr<NпјҢжҲ‘жғід»ҺдёӯжҸҗеҸ–rдёӘж–№зЁӢпјҢд»ҚеҢ…еҗ«е…ЁйғЁдҝЎжҒҜгҖӮ - flonk

你可以在这里查看：http://docs.sympy.org/dev/modules/solvers/solvers.html - YXD

谢谢E先生，但如果可能的话，我想避免将其转换为符号对象。最好有一些明确针对数组（浮点数）的东西，在最好的情况下，对于整数数组变得精确。 - flonk

我自己写了一个 这里 - ricardo

Answer 1

我终于发现，可以使用LU分解来完成。这里的U矩阵代表线性系统的简化形式。

from numpy import array
from scipy.linalg import lu

a = array([[2.,4.,4.,4.],[1.,2.,3.,3.],[1.,2.,2.,2.],[1.,4.,3.,4.]])

pl, u = lu(a, permute_l=True)

那么u读取

array([[ 2.,  4.,  4.,  4.],
       [ 0.,  2.,  1.,  2.],
       [ 0.,  0.,  1.,  1.],
       [ 0.,  0.,  0.,  0.]])

根据系统可解性，矩阵具有上三角形或梯形结构。在上述情况中，会出现一行零，因为该矩阵仅具有秩为3。

- flonk

---

请看这个例子：http://poquitopicante.blogspot.com/2014/04/gaussian-elimination-using-lu.html - Mark Mikofski

@alvas 矩阵的LU分解是将其分解为一个下三角矩阵和一个上三角矩阵的乘积，参见例如http://en.wikipedia.org/wiki/LU_decomposition。 - flonk

在Lapack中，U因子不需要处于行梯阵形式。请参见https://math.stackexchange.com/questions/891334/is-the-u-factor-in-lu-decomposition-for-rectangular-matrices-always-in-row-echel - Daniel Arteaga

Answer 2

如果您需要删除重复或冗余的等式，可能值得检查的一个函数是_remove_redundancy。

import numpy as np
import scipy.optimize

a = np.array([[1.,1.,1.,1.],
              [0.,0.,0.,1.],
              [0.,0.,0.,2.],
              [0.,0.,0.,3.]])
print(scipy.optimize._remove_redundancy._remove_redundancy(a, np.zeros_like(a[:, 0]))[0])

这将会得到：

[[1. 1. 1. 1.]
 [0. 0. 0. 3.]]

作为对@flonk答案的补充说明，使用LU分解可能不能始终得到所需的降阶矩阵。例如：

import numpy as np
import scipy.linalg

a = np.array([[1.,1.,1.,1.],
              [0.,0.,0.,1.],
              [0.,0.,0.,2.],
              [0.,0.,0.,3.]])

_,_, u = scipy.linalg.lu(a)
print(u)

生成相同的矩阵：

[[1. 1. 1. 1.]
 [0. 0. 0. 1.]
 [0. 0. 0. 2.]
 [0. 0. 0. 3.]]

即使最后3行线性相关。

- Bernardo Costa

Answer 3

您可以使用Python符号数学库sympy。

import sympy as sp

m = sp.Matrix([[1,2,1],
           [-2,-3,1],
           [3,5,0]])

m_rref, pivots = m.rref() # Compute reduced row echelon form (rref).

print(m_rref, pivots)

此命令将输出矩阵的简化阶梯形式，以及主元列的列表。

Matrix([[1, 0, -5],
        [0, 1,  3],
        [0, 0,  0]])

(0, 1)

- rsant