Reduced row echelon form

rref.py

import numpy as np


def reduced_row_echelon_form(matrix: np.ndarray, atol=1e-8):
    """
    return reduced row echelon form of a given matrix

    Parameters
    ----------
    matrix: (nr, nl)
    atol: absolute tolerance for float comparison

    Returns
    -------
    rref: (nr, nl)
        reduced row echelon form
    transformation: (nr, nr)
        rref = np.dot(transformation, matrix)
    """
    assert matrix.ndim == 2
    nr, nl = matrix.shape

    rref = matrix.copy()
    transformation = np.eye(nr)

    def swap_i_j(rref, transformation, i, j):
        # swap axis-i and axis-j
        assert i != j
        rref[[i, j]] = rref[[j, i]]

        mat = np.eye(nr)
        mat[i, i] = mat[j, j] = 0
        mat[i, j] = mat[j, i] = 1
        transformation = np.dot(mat, transformation)
        return rref, transformation

    def add_i_jk(rref, transformation, i, j, k):
        # add k-times axis-j to axis-i
        assert i != j
        rref[i] += k * rref[j]

        mat = np.eye(nr)
        mat[i, j] = k
        transformation = np.dot(mat, transformation)
        return rref, transformation

    def multiply_i(rref, transformation, i, k):
        # multiply axis-i by k
        rref[i] *= k

        mat = np.eye(nr)
        mat[i, i] = k
        transformation = np.dot(mat, transformation)
        return rref, transformation

    def iszero(x: float) -> bool:
        return np.isclose(x, 0., atol=atol)

    pivot = 0
    for row in range(nr):
        if pivot >= nl:
            break

        # determine pivot
        found_nonzero = False
        for i in range(row, nr):
            if not iszero(rref[i, pivot]):
                if i != row:
                    rref, transformation = swap_i_j(rref, transformation, row, i)
                found_nonzero = True
                break

        if found_nonzero:
            # here, rref[row, pivot] != 0
            if not np.isclose(rref[row, pivot], 1, atol=atol):
                rref, transformation = multiply_i(rref, transformation, row, 1. / rref[row, pivot])

            for i in range(nr):
                if (i != row) and not iszero(rref[i, pivot]):
                    rref, transformation = add_i_jk(rref, transformation, i, row, -rref[i, pivot])
            pivot += 1

    return rref, transformation


def is_row_echelon_form(matrix, rref, transformation, atol=1e-8) -> bool:
    assert matrix.ndim == 2
    nr, nl = matrix.shape
    assert transformation.shape == (nr, nr)
    return np.allclose(rref, np.dot(transformation, matrix), atol=atol)


if __name__ == '__main__':
    matrix = np.array([
        [0, 1, 1],
        [2, 2, 4],
    ], dtype=np.float)

    transformation = np.array([
        [-1, 0.5],
        [1, 0],
    ])
    rref = np.array([
        [1, 0, 1],
        [0, 1, 1],
    ], dtype=np.float)
    assert is_row_echelon_form(matrix, rref, transformation)

    rref_expect, transformation_expect = reduced_row_echelon_form(matrix)
    assert np.allclose(rref_expect, rref)
    assert np.allclose(transformation_expect, transformation)

    matrix = np.array([
        [ 1.,  0.,  0.,  0.],
        [ 0.,  0.,  0.,  1.],
        [ 0.,  1.,  0.,  0.],
        [-3., -1.,  1.,  1.],
        [-2., -1.,  2.,  0.],
        [ 0.,  0., -3.,  1.]
    ])
    rref, transformation = reduced_row_echelon_form(matrix)
    assert is_row_echelon_form(matrix, rref, transformation)