Recursive Gauss Elimination algorithm

Question

I'm studying functional programming concepts in Python 3, so I wrote this Gauss Elimination algorithm with tail-recursion. The algorithm passes all tests I found in a text-book except for one. I had to write simple Matrix and Vector classes along with some accessory functions, because the web-platform I use doesn't have NumPy (the classes are provided at the bottom).

NB! Based on the first answer I received I should stress that there is a cutoff for arithmetics precision (you can find DIGITS in the block at the bottom)

The algorithm

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

Tests

# 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)))

Let's see how it transforms the matrix step by step:

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

I get a consistent system with infinite number of solutions. Everything seems right to me, yet I'm supposed to get this inconsistent system (I checked it with NumPy and R):

((1, 0, -3/23, -1/23, 0),
 (0, 1, -11/23, -19/23, 0),
 (0, 0, 0, 0, 1)) 

I hope I'm missing something obvious and simple. I hope the code is documented well enough to be readable. Thank you in advance.

Classes, accessory functions and constants:

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))

Answer 1

I see some very large and very small numbers around the last steps; it seems you're victim of rounding errors, producing not 0 but very small numbers (and their very large inverses).

You would need to set a (relative) cutoff below which a number is deemed zero (2e-15 seems a good one). Find the epsilon for your system: the minimum difference between two floats that can be properly represented and are of order 1.

Do check methods like normalise_row(), where you do

coeff = 1 / (self[i][j]) if self[i][j] else 1

and a few other places where I don't see a use of round().

Another possibility might be to use the Python decimal module. I have not much experience with it, but perhaps it provides the necessary precision for your calculations.

Answer 2

I only had a quick look at the algorithm and it seems to be doing what you would expect of it. The real problem is in the "mechanics" of how you would detect an inconsistent solution. If you look at the 6th step of your algorithm you get

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

Note that the 3rd and 4th terms on the third line are very tiny compared with the other numbers. In fact they are just round off errors. If you did the subtraction with exact numbers then the 2nd column would give you 46/3 - 2*23/2 = 0 (as you get), but the third column would give you -22/3 - 2*(-11/3) = 0 as you don't get. Similarly for the fourth column. Then in the next stage you compound the error by "scaling" these zeros up to 1 and 2 and all will be in error from there onward.

You have two choices - either do the calculations with exact numbers (I don't know if python can do this but, for example, in Scheme there are exact fractions that would not produce this problem - maybe you would have to roll your own...) - or alternately, use a sufficiently high precision arithmetic and arbitrarily set anything smaller than some chosen value (best based on the magnitude of the numbers in the original matrix) equal to zero. This would fail if the matrix really did have a solution with numbers of wildly varying magnitude, but in 'real-world' problems it should weed out the problems you are having.