递归高斯消除算法
Recursive Gauss Elimination algorithm
我正在研究 Python 3 中的函数式编程概念，所以我用尾递归编写了这个高斯消元算法。该算法通过了我在教科书中找到的所有测试，除了一个。我不得不编写简单的 Matrix 和 Vector 类 以及一些辅助函数，因为我使用的网络平台没有 NumPy（类 在底部提供）。
注意！根据我收到的第一个答案，我应该强调算术精度有一个截止点（你可以在底部的块中找到 DIGITS）
算法
def istrue(iterable, key=bool):
"""
:type iterable: collections.Iterable
:type key: (object) -> bool
:rtype: tuple[bool]
"""
try:
return tuple(map(int, map(key, iterable)))
except TypeError:
return bool(iterable),
def where(iterable, key=bool):
"""
:type iterable: collections.Iterable
:type key: (object) -> bool
:rtype: tuple[int]
"""
return tuple(i for i, elem in enumerate(iterable) if key(elem))
def next_true(iterable, i, key=bool):
"""
Returns position of a True element next to the i'th element
(iterable[j]: j>i, key(iterable[j]) -> True) if it exists.
:type iterable: collections.Iterable
:type i: int
:rtype: int | None
"""
true_mask = istrue(iterable, key)
try:
return where(enumerate(true_mask),
key=lambda ind_val: ind_val[0] > i and ind_val[1])[0]
except IndexError:
# No elements satisfy the condition
return None
def row_echelon(matrix):
"""
Transforms matrix into the row echelon form
:type matrix: Matrix
:return: upper-triangular matrix
:rtype: Matrix
"""
@optimize_tail_recursion
def operator(m, var):
"""
:type m: Matrix
:type var: int
:rtype: Matrix
"""
# if we have already processed all variables we could
if var > m.ncol - 1 or var > m.nrow - 1:
return m
# if matrix[var][var] is zero and there exist i such that
# matrix[i][var] > 0, i > j
elif not m[var][var] and sum(istrue(m.col(var))[var:]):
i = next_true(istrue(m.col(var)), var)
return operator(m.permute(var, i), var)
# if |{matrix[i][var], 0 <= i < nrow(matrix)}| > 1, matrix[var][var]!=0
elif m[var][var] and sum(istrue(m.col(var))) - 1:
i = tuple(ind for ind in where(m.col(var)) if ind != var)[0]
coeff = - m[i][var] / m[var][var]
return operator(m.add_rows(var, i, coeff), var)
# if matrix[var][var] is zero and there is no i such that
# matrix[i][var] > 0, i > j
# or all possible elements were eliminated
return operator(m.normalise_row(var, var), var+1)
return operator(matrix, 0)
def solve_linear_system(augmented_matrix):
"""
:type augmented_matrix: Matrix
:return: str | tuple[float]
"""
row_echelon_m = row_echelon(augmented_matrix)
print(row_echelon_m)
left_side, right_side = (row_echelon_m[:, :row_echelon_m.ncol-1],
row_echelon_m.col(-1))
nontrivial = istrue(left_side, key=lambda row: bool(sum(row)))
rank = sum(nontrivial)
# if there exists an equation with trivial coefficients and nontrivial
# right-hand side
if sum(not sum(l_side) and r_side for l_side, r_side
in zip(left_side, right_side)):
return NO
# rank(matrix) < number of variables
elif rank < left_side.ncol:
return INF
return right_side
测试
# must be (10, 5, -20.0) - passed
test1 = Matrix(((0.12, 0.18, -0.17, 5.5), (0.06, 0.09, 0.15, -1.95), (0.22, -0.1, 0.06, 0.5)))
# infinite number of solution - passed
test2 = Matrix((((1, 2, -1, 3, 7), (2, 4, -2, 6, 14), (1, -1, 3, 1, -1)))
# no solutions - passed
test3 = Matrix(((2, -1, 3, 1), (2, -1, -1, -2), (4, -2, 6, 0), (6, 8, -7, 2)))
# infinite number of solution - passed
test4 = Matrix(((3, -2, 1, 0), (5, -14, 15, 0), (1, 2, -3, 0)))
# infinite number of solution - passed
test5 = Matrix((((2, 3, -1, 1, 0), (2, 7, -3, 0, 1), (0, 4, -2, -1, 1), (2, -1, 1, 2, -1), (4, 10, -4, 1, 1)))
# no solutions - failed. My algorithm returns Inf
test6 = Matrix(((3, -5, 2, 4, 2), (7, -4, 1, 3, 5), (5, 7, -4, -6, 3)))
来看看它是如何一步步变换矩阵的：
solve_linear_system(test6) # ->
# The integer on the right side denoted the variable that is being eliminated at the moment using the corresponding row.
((3.0, -5.0, 2.0, 4.0, 2.0),
(7.0, -4.0, 1.0, 3.0, 5.0),
(5.0, 7.0, -4.0, -6.0, 3.0)) 0
((3.0, -5.0, 2.0, 4.0, 2.0),
(0.0, 7.666666666666668, -3.666666666666667, -6.333333333333334, 0.33333333333333304),
(5.0, 7.0, -4.0, -6.0, 3.0)) 0
((3.0, -5.0, 2.0, 4.0, 2.0),
(0.0, 7.666666666666668, -3.666666666666667, -6.333333333333334, 0.33333333333333304),
(0.0, 15.333333333333334, -7.333333333333334, -12.666666666666668, -0.3333333333333335)) 0
((1.0, -1.6666666666666665, 0.6666666666666666, 1.3333333333333333, 0.6666666666666666),
(0.0, 7.666666666666668, -3.666666666666667, -6.333333333333334, 0.33333333333333304),
(0.0, 15.333333333333334, -7.333333333333334, -12.666666666666668, -0.3333333333333335)) 1
((1.0, 0.0, -0.13043478260869557, -0.043478260869564966, 0.7391304347826085),
(0.0, 7.666666666666668, -3.666666666666667, -6.333333333333334, 0.33333333333333304),
(0.0, 15.333333333333334, -7.333333333333334, -12.666666666666668, -0.3333333333333335)) 1
((1.0, 0.0, -0.13043478260869557, -0.043478260869564966, 0.7391304347826085),
(0.0, 7.666666666666668, -3.666666666666667, -6.333333333333334, 0.33333333333333304),
(0.0, 0.0, -8.88e-16, -1.776e-15, -0.9999999999999994)) 1
((1.0, 0.0, -0.13043478260869557, -0.043478260869564966, 0.7391304347826085),
(0.0, 0.9999999999999999, -0.4782608695652173, -0.826086956521739, 0.04347826086956517),
(0.0, 0.0, -8.88e-16, -1.776e-15, -0.9999999999999994)) 2
((1.0, 0.0, 0.0, 0.21739130434782616, 146886016451234.4),
(0.0, 0.9999999999999999, -0.4782608695652173, -0.826086956521739, 0.04347826086956517),
(0.0, 0.0, -8.88e-16, -1.776e-15, -0.9999999999999994)) 2
((1.0, 0.0, 0.0, 0.21739130434782616, 146886016451234.4),
(0.0, 0.9999999999999999, 0.0, 0.13043478260869557, 538582060321190.4),
(0.0, 0.0, -8.88e-16, -1.776e-15, -0.9999999999999994)) 2
((1.0, 0.0, 0.0, 0.21739130434782616, 146886016451234.4),
(0.0, 0.9999999999999999, 0.0, 0.13043478260869557, 538582060321190.4),
(-0.0, -0.0, 0.9999999999999999, 1.9999999999999998, 1126126126126125.2)) 3
我得到了一个具有无限数量解决方案的一致系统。一切对我来说似乎都是正确的，但我应该得到这个不一致的系统（我用 NumPy 和 R 检查了它）：
((1, 0, -3/23, -1/23, 0),
(0, 1, -11/23, -19/23, 0),
(0, 0, 0, 0, 1))
我希望我遗漏了一些明显而简单的东西。我希望代码有足够的文档记录以便于阅读。先感谢您。
类，附属函数和常量：
NO = "NO"
YES = "YES"
INF = "INF"
DIGITS = 16
def typemap(iterable, types):
try:
return all(isinstance(elem, _type)
for elem, _type in zip(iterable, types))
except TypeError:
return False
class TailRecursionException(BaseException):
def __init__(self, args, kwargs):
self.args = args
self.kwargs = kwargs
def optimize_tail_recursion(g):
def func(*args, **kwargs):
f = sys._getframe()
if f.f_back and f.f_back.f_back and f.f_back.f_back.f_code == f.f_code:
raise TailRecursionException(args, kwargs)
else:
while 1:
try:
return g(*args, **kwargs)
except TailRecursionException as e:
args = e.args
kwargs = e.kwargs
func.__doc__ = g.__doc__
return func
class Vector(tuple):
def __add__(self, other):
if not isinstance(other, tuple):
raise TypeError
return Vector(round(a + b, DIGITS) for a, b in zip(self, other))
def __radd__(self, other):
return self.__add__(other)
def __sub__(self, other):
return self.__add__(-1 * other)
def __rsub__(self, other):
return self.__sub__(other)
def __mul__(self, other):
if not isinstance(other, (int, float)):
raise TypeError
return Vector(round(elem * round(other, DIGITS), DIGITS) for elem in self)
def __rmul__(self, other):
return self.__mul__(other)
def extend(self, item):
return Vector(super().__add__((item, )))
def concat(self, item):
return Vector(super().__add__(item))
class Matrix:
"""
:type _matrix: tuple[Vector]
"""
def __init__(self, matrix):
"""
:type matrix: list[list] | tuple[tuple]
:return:
"""
if not matrix:
raise ValueError("Empty matrices are not supported")
self._matrix = tuple(map(Vector, map(lambda row: map(float, row),
matrix)))
def __iter__(self):
return iter(self._matrix)
def __repr__(self):
return "({})".format(",\n ".join(map(repr, self)))
def __getitem__(self, item):
"""
Returns a Vector if item is int or slice; returns a Matrix if item is
tuple of slices
:param item:
:rtype: Matrix | Vector
"""
if isinstance(item, (int, slice)):
return self._matrix[item]
elif typemap(item, (slice, slice)):
row_slice, col_slice = item
return Matrix(tuple(map(op.itemgetter(col_slice), self[row_slice])))
def __mul__(self, coeff):
if not isinstance(coeff, (int, float)):
raise TypeError
return Matrix(tuple(vector * coeff for vector in self))
def __rmul__(self, coeff):
return self.__mul__(coeff)
@property
def nrow(self):
return len(self._matrix)
@property
def ncol(self):
return len(self.row(0))
def row(self, i):
return self[i]
def col(self, j):
"""
:type j: int
:return: tuple[tuple]
"""
return Vector(self[i][j] for i in range(self.nrow))
def _replace_row(self, i, replacement):
new_matrix = tuple(self.row(_i) if _i != i else replacement
for _i in range(self.nrow))
return Matrix(new_matrix)
def permute(self, a, b):
"""
Exchange rows a and b
:type a: int
:type b: int
:rtype: Matrix
"""
return self._replace_row(b, self.row(a))._replace_row(a, self.row(b))
def multiply_row(self, i, coeff):
return self._replace_row(i, self.row(i) * coeff)
def normalise_row(self, i, j):
coeff = 1 / (self[i][j]) if self[i][j] else 1
return self.multiply_row(i, coeff)
def add_rows(self, a, b, coeff=1):
"""
:return: matrix': matrix'[b] = coef * matrix[a] + matrix[b]
:rtype: Matrix
"""
return self._replace_row(b, coeff * self.row(a) + self.row(b))
我在最后几步看到一些非常大和非常小的数字；看来您是舍入错误的受害者，产生的不是 0，而是非常小的数字（及其非常大的倒数）。
您需要设置一个（相对）截止值，低于该值的数字被视为零（2e-15 似乎是一个不错的值）。找到适合您系统的 epsilon：可以正确表示且阶数为 1 的两个浮点数之间的最小差异。
检查 normalise_row() 之类的方法，
coeff = 1 / (self[i][j]) if self[i][j] else 1
以及其他一些我没有看到 round().
用法的地方
另一种可能是使用 Python decimal 模块。我对此经验不多，但也许它为您的计算提供了必要的精度。
我只是快速浏览了一下算法，它似乎在做您期望的事情。真正的问题在于 "mechanics" 如何检测不一致的解决方案。如果您查看算法的第 6 步，您会得到
((1.0, 0.0, -0.13043478260869557, -0.043478260869564966, 0.7391304347826085),
(0.0, 7.666666666666668, -3.666666666666667, -6.333333333333334, 0.33333333333333304),
(0.0, 0.0, -8.88e-16, -1.776e-15, -0.9999999999999994))
请注意，与其他数字相比，第三行的第 3 项和第 4 项非常小。事实上，它们只是四舍五入的错误。如果你用精确的数字做减法，那么第二列会给你 46/3 - 2*23/2 = 0 （正如你得到的那样），但第三列会给你 -22/3 - 2*(-11/ 3) = 0 因为你没有得到。第四列也类似。然后在下一阶段，你通过 "scaling" 这些零到 1 和 2 来复合错误，并且从那里开始所有错误。
你有两个选择 - 要么用精确的数字进行计算（我不知道 python 是否可以这样做，但是，例如，在 Scheme 中有精确的分数不会产生这个问题 -也许你必须自己动手......） - 或者，使用足够高的精度算术并任意将小于某个选定值（最好基于原始矩阵中数字的大小）的任何值设置为零。如果矩阵确实有一个数量级变化很大的解决方案，这将失败，但在 'real-world' 问题中，它应该消除你遇到的问题。
我正在研究 Python 3 中的函数式编程概念，所以我用尾递归编写了这个高斯消元算法。该算法通过了我在教科书中找到的所有测试，除了一个。我不得不编写简单的 Matrix 和 Vector 类 以及一些辅助函数，因为我使用的网络平台没有 NumPy（类 在底部提供）。
注意！根据我收到的第一个答案，我应该强调算术精度有一个截止点（你可以在底部的块中找到 DIGITS）
算法
def istrue(iterable, key=bool):
"""
:type iterable: collections.Iterable
:type key: (object) -> bool
:rtype: tuple[bool]
"""
try:
return tuple(map(int, map(key, iterable)))
except TypeError:
return bool(iterable),
def where(iterable, key=bool):
"""
:type iterable: collections.Iterable
:type key: (object) -> bool
:rtype: tuple[int]
"""
return tuple(i for i, elem in enumerate(iterable) if key(elem))
def next_true(iterable, i, key=bool):
"""
Returns position of a True element next to the i'th element
(iterable[j]: j>i, key(iterable[j]) -> True) if it exists.
:type iterable: collections.Iterable
:type i: int
:rtype: int | None
"""
true_mask = istrue(iterable, key)
try:
return where(enumerate(true_mask),
key=lambda ind_val: ind_val[0] > i and ind_val[1])[0]
except IndexError:
# No elements satisfy the condition
return None
def row_echelon(matrix):
"""
Transforms matrix into the row echelon form
:type matrix: Matrix
:return: upper-triangular matrix
:rtype: Matrix
"""
@optimize_tail_recursion
def operator(m, var):
"""
:type m: Matrix
:type var: int
:rtype: Matrix
"""
# if we have already processed all variables we could
if var > m.ncol - 1 or var > m.nrow - 1:
return m
# if matrix[var][var] is zero and there exist i such that
# matrix[i][var] > 0, i > j
elif not m[var][var] and sum(istrue(m.col(var))[var:]):
i = next_true(istrue(m.col(var)), var)
return operator(m.permute(var, i), var)
# if |{matrix[i][var], 0 <= i < nrow(matrix)}| > 1, matrix[var][var]!=0
elif m[var][var] and sum(istrue(m.col(var))) - 1:
i = tuple(ind for ind in where(m.col(var)) if ind != var)[0]
coeff = - m[i][var] / m[var][var]
return operator(m.add_rows(var, i, coeff), var)
# if matrix[var][var] is zero and there is no i such that
# matrix[i][var] > 0, i > j
# or all possible elements were eliminated
return operator(m.normalise_row(var, var), var+1)
return operator(matrix, 0)
def solve_linear_system(augmented_matrix):
"""
:type augmented_matrix: Matrix
:return: str | tuple[float]
"""
row_echelon_m = row_echelon(augmented_matrix)
print(row_echelon_m)
left_side, right_side = (row_echelon_m[:, :row_echelon_m.ncol-1],
row_echelon_m.col(-1))
nontrivial = istrue(left_side, key=lambda row: bool(sum(row)))
rank = sum(nontrivial)
# if there exists an equation with trivial coefficients and nontrivial
# right-hand side
if sum(not sum(l_side) and r_side for l_side, r_side
in zip(left_side, right_side)):
return NO
# rank(matrix) < number of variables
elif rank < left_side.ncol:
return INF
return right_side
测试
# must be (10, 5, -20.0) - passed
test1 = Matrix(((0.12, 0.18, -0.17, 5.5), (0.06, 0.09, 0.15, -1.95), (0.22, -0.1, 0.06, 0.5)))
# infinite number of solution - passed
test2 = Matrix((((1, 2, -1, 3, 7), (2, 4, -2, 6, 14), (1, -1, 3, 1, -1)))
# no solutions - passed
test3 = Matrix(((2, -1, 3, 1), (2, -1, -1, -2), (4, -2, 6, 0), (6, 8, -7, 2)))
# infinite number of solution - passed
test4 = Matrix(((3, -2, 1, 0), (5, -14, 15, 0), (1, 2, -3, 0)))
# infinite number of solution - passed
test5 = Matrix((((2, 3, -1, 1, 0), (2, 7, -3, 0, 1), (0, 4, -2, -1, 1), (2, -1, 1, 2, -1), (4, 10, -4, 1, 1)))
# no solutions - failed. My algorithm returns Inf
test6 = Matrix(((3, -5, 2, 4, 2), (7, -4, 1, 3, 5), (5, 7, -4, -6, 3)))
来看看它是如何一步步变换矩阵的：
solve_linear_system(test6) # ->
# The integer on the right side denoted the variable that is being eliminated at the moment using the corresponding row.
((3.0, -5.0, 2.0, 4.0, 2.0),
(7.0, -4.0, 1.0, 3.0, 5.0),
(5.0, 7.0, -4.0, -6.0, 3.0)) 0
((3.0, -5.0, 2.0, 4.0, 2.0),
(0.0, 7.666666666666668, -3.666666666666667, -6.333333333333334, 0.33333333333333304),
(5.0, 7.0, -4.0, -6.0, 3.0)) 0
((3.0, -5.0, 2.0, 4.0, 2.0),
(0.0, 7.666666666666668, -3.666666666666667, -6.333333333333334, 0.33333333333333304),
(0.0, 15.333333333333334, -7.333333333333334, -12.666666666666668, -0.3333333333333335)) 0
((1.0, -1.6666666666666665, 0.6666666666666666, 1.3333333333333333, 0.6666666666666666),
(0.0, 7.666666666666668, -3.666666666666667, -6.333333333333334, 0.33333333333333304),
(0.0, 15.333333333333334, -7.333333333333334, -12.666666666666668, -0.3333333333333335)) 1
((1.0, 0.0, -0.13043478260869557, -0.043478260869564966, 0.7391304347826085),
(0.0, 7.666666666666668, -3.666666666666667, -6.333333333333334, 0.33333333333333304),
(0.0, 15.333333333333334, -7.333333333333334, -12.666666666666668, -0.3333333333333335)) 1
((1.0, 0.0, -0.13043478260869557, -0.043478260869564966, 0.7391304347826085),
(0.0, 7.666666666666668, -3.666666666666667, -6.333333333333334, 0.33333333333333304),
(0.0, 0.0, -8.88e-16, -1.776e-15, -0.9999999999999994)) 1
((1.0, 0.0, -0.13043478260869557, -0.043478260869564966, 0.7391304347826085),
(0.0, 0.9999999999999999, -0.4782608695652173, -0.826086956521739, 0.04347826086956517),
(0.0, 0.0, -8.88e-16, -1.776e-15, -0.9999999999999994)) 2
((1.0, 0.0, 0.0, 0.21739130434782616, 146886016451234.4),
(0.0, 0.9999999999999999, -0.4782608695652173, -0.826086956521739, 0.04347826086956517),
(0.0, 0.0, -8.88e-16, -1.776e-15, -0.9999999999999994)) 2
((1.0, 0.0, 0.0, 0.21739130434782616, 146886016451234.4),
(0.0, 0.9999999999999999, 0.0, 0.13043478260869557, 538582060321190.4),
(0.0, 0.0, -8.88e-16, -1.776e-15, -0.9999999999999994)) 2
((1.0, 0.0, 0.0, 0.21739130434782616, 146886016451234.4),
(0.0, 0.9999999999999999, 0.0, 0.13043478260869557, 538582060321190.4),
(-0.0, -0.0, 0.9999999999999999, 1.9999999999999998, 1126126126126125.2)) 3
我得到了一个具有无限数量解决方案的一致系统。一切对我来说似乎都是正确的，但我应该得到这个不一致的系统（我用 NumPy 和 R 检查了它）：
((1, 0, -3/23, -1/23, 0),
(0, 1, -11/23, -19/23, 0),
(0, 0, 0, 0, 1))
我希望我遗漏了一些明显而简单的东西。我希望代码有足够的文档记录以便于阅读。先感谢您。
类，附属函数和常量：
NO = "NO"
YES = "YES"
INF = "INF"
DIGITS = 16
def typemap(iterable, types):
try:
return all(isinstance(elem, _type)
for elem, _type in zip(iterable, types))
except TypeError:
return False
class TailRecursionException(BaseException):
def __init__(self, args, kwargs):
self.args = args
self.kwargs = kwargs
def optimize_tail_recursion(g):
def func(*args, **kwargs):
f = sys._getframe()
if f.f_back and f.f_back.f_back and f.f_back.f_back.f_code == f.f_code:
raise TailRecursionException(args, kwargs)
else:
while 1:
try:
return g(*args, **kwargs)
except TailRecursionException as e:
args = e.args
kwargs = e.kwargs
func.__doc__ = g.__doc__
return func
class Vector(tuple):
def __add__(self, other):
if not isinstance(other, tuple):
raise TypeError
return Vector(round(a + b, DIGITS) for a, b in zip(self, other))
def __radd__(self, other):
return self.__add__(other)
def __sub__(self, other):
return self.__add__(-1 * other)
def __rsub__(self, other):
return self.__sub__(other)
def __mul__(self, other):
if not isinstance(other, (int, float)):
raise TypeError
return Vector(round(elem * round(other, DIGITS), DIGITS) for elem in self)
def __rmul__(self, other):
return self.__mul__(other)
def extend(self, item):
return Vector(super().__add__((item, )))
def concat(self, item):
return Vector(super().__add__(item))
class Matrix:
"""
:type _matrix: tuple[Vector]
"""
def __init__(self, matrix):
"""
:type matrix: list[list] | tuple[tuple]
:return:
"""
if not matrix:
raise ValueError("Empty matrices are not supported")
self._matrix = tuple(map(Vector, map(lambda row: map(float, row),
matrix)))
def __iter__(self):
return iter(self._matrix)
def __repr__(self):
return "({})".format(",\n ".join(map(repr, self)))
def __getitem__(self, item):
"""
Returns a Vector if item is int or slice; returns a Matrix if item is
tuple of slices
:param item:
:rtype: Matrix | Vector
"""
if isinstance(item, (int, slice)):
return self._matrix[item]
elif typemap(item, (slice, slice)):
row_slice, col_slice = item
return Matrix(tuple(map(op.itemgetter(col_slice), self[row_slice])))
def __mul__(self, coeff):
if not isinstance(coeff, (int, float)):
raise TypeError
return Matrix(tuple(vector * coeff for vector in self))
def __rmul__(self, coeff):
return self.__mul__(coeff)
@property
def nrow(self):
return len(self._matrix)
@property
def ncol(self):
return len(self.row(0))
def row(self, i):
return self[i]
def col(self, j):
"""
:type j: int
:return: tuple[tuple]
"""
return Vector(self[i][j] for i in range(self.nrow))
def _replace_row(self, i, replacement):
new_matrix = tuple(self.row(_i) if _i != i else replacement
for _i in range(self.nrow))
return Matrix(new_matrix)
def permute(self, a, b):
"""
Exchange rows a and b
:type a: int
:type b: int
:rtype: Matrix
"""
return self._replace_row(b, self.row(a))._replace_row(a, self.row(b))
def multiply_row(self, i, coeff):
return self._replace_row(i, self.row(i) * coeff)
def normalise_row(self, i, j):
coeff = 1 / (self[i][j]) if self[i][j] else 1
return self.multiply_row(i, coeff)
def add_rows(self, a, b, coeff=1):
"""
:return: matrix': matrix'[b] = coef * matrix[a] + matrix[b]
:rtype: Matrix
"""
return self._replace_row(b, coeff * self.row(a) + self.row(b))
我在最后几步看到一些非常大和非常小的数字；看来您是舍入错误的受害者，产生的不是 0，而是非常小的数字（及其非常大的倒数）。
您需要设置一个（相对）截止值，低于该值的数字被视为零（2e-15 似乎是一个不错的值）。找到适合您系统的 epsilon：可以正确表示且阶数为 1 的两个浮点数之间的最小差异。
检查 normalise_row() 之类的方法，
coeff = 1 / (self[i][j]) if self[i][j] else 1
以及其他一些我没有看到 round().
另一种可能是使用 Python decimal 模块。我对此经验不多，但也许它为您的计算提供了必要的精度。
我只是快速浏览了一下算法，它似乎在做您期望的事情。真正的问题在于 "mechanics" 如何检测不一致的解决方案。如果您查看算法的第 6 步，您会得到
((1.0, 0.0, -0.13043478260869557, -0.043478260869564966, 0.7391304347826085),
(0.0, 7.666666666666668, -3.666666666666667, -6.333333333333334, 0.33333333333333304),
(0.0, 0.0, -8.88e-16, -1.776e-15, -0.9999999999999994))
请注意，与其他数字相比，第三行的第 3 项和第 4 项非常小。事实上，它们只是四舍五入的错误。如果你用精确的数字做减法，那么第二列会给你 46/3 - 2*23/2 = 0 （正如你得到的那样），但第三列会给你 -22/3 - 2*(-11/ 3) = 0 因为你没有得到。第四列也类似。然后在下一阶段，你通过 "scaling" 这些零到 1 和 2 来复合错误，并且从那里开始所有错误。
你有两个选择 - 要么用精确的数字进行计算（我不知道 python 是否可以这样做，但是，例如，在 Scheme 中有精确的分数不会产生这个问题 -也许你必须自己动手......） - 或者，使用足够高的精度算术并任意将小于某个选定值（最好基于原始矩阵中数字的大小）的任何值设置为零。如果矩阵确实有一个数量级变化很大的解决方案，这将失败，但在 'real-world' 问题中，它应该消除你遇到的问题。