Using Python and NumPy for Computations in Linear Algebra¶

import numpy as np
import numpy.linalg as la

A = np.random.randint(-10, 10, size=(2,3))
B = np.random.randint(-10, 10, size=(2,3))

A = np.array([2, 8, 4, 2, 5, 1, 4, 10, -1]).reshape(3, 3)
B = np.array([1, 2, 0, 1, 1, 1]).reshape(3, 2)

def clear_col(M, i, j, precision):
    """
    Given matrix M and position (i, j), if M[i, j] != 0, then make it 1,
    and make all entries above and below it equal to 0.
    Since we are dearling with floats, add a precision for comparison.
    """
    if M[i, j] == 0:
        raise ValueError("Element equal to zero!")

    # copy the matrix, as to non change the original!
    N = M.copy()

    # make entry 1
    N[i] *= 1/N[i, j]

    # clear
    for row in range(N.shape[0]):
        if row == i or np.abs(N[row, j]) < precision:
            continue
        N[row] += -N[row, j] * N[i]
    return N


def rref(M, precision=1e-9, type=np.float64, lastcol=None):
    """
    Given matrix M and a last column lastcol, perform row operations
    to put only up to the lastcol part in reduced row echelon form.
    Since we are performing division, we need to convert to floats.  (Defaut: np.float64)
    Since we are dearling with floats, add a precision for comparison.  (Default: 9 digits)
    """
    # make sure the entries are floats (also makes unnecessary to use copy)
    N = M.astype(type)

    # Make the default last column the real last one
    if lastcol is None:
        lastcol = N.shape[1]

    # main part
    cur_row = 0
    n_rows = N.shape[0]
    for j in range(lastcol):
        for i in range(cur_row, n_rows):
            if np.abs(N[i, j]) > precision:
                if i != cur_row:
                    # N[i], N[cur_row] = N[cur_row], N[i]
                    N[[i, cur_row]] = N[[cur_row, i]]
                N = clear_col(N, cur_row, j, precision)
                cur_row += 1
    return N

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

def remove_zero_rows(A, precision=1e-9):
    """
    Remove rows of zeros.
    """
    # produce a "mask": True for non-zero entries, False for zero entries
    TF = np.abs(A) > precision

    # return rows where we have a True entry
    return A[np.any(TF, axis=1)]


def find_leaing_ones(E, precision=1e-9):
    """
    Given matrix E in reduced row echelon form, returns a list with the positions of leading ones.
    """
    # produce a "mask": True for non-zero entries, False for zero entries
    TF = np.abs(E) > precision

    # drop rows with only False
    TF = TF[np.any(TF, axis=1)]

    res = []
    for i, row in enumerate(TF):
        # argmax gives the index/argument for first True (leading 1 position) of the row
        res.append((i, np.argmax(row)))
    return res


def rank(A, precision=1e-9):
    """
    Rank of a given matrix A.
    """
    return len(find_leaing_ones(rref(A, precision=precision), precision=precision))


def nullity(A, precision=1e-9):
    """
    Nullitty of given matrix A.
    """
    return A.shape[1] - rank(A, precision=precision)


def row_space(A, precision=1e-9):
    """
    Given a matrix A, reutrns the basis of the row space of A, as rows of a matrix.
    """
    return remove_zero_rows(rref(A, precision=precision), precision=precision)


def column_space(A, precision=1e-9):
    """
    Given a matrix A, reutrns the basis of the column space of A, as *rows* of a matrix.
    """
    return row_space(A.T, precision=precision)


def nullspace(A, precision=1e-9):
    """
    Given A, returns a basis for the nullspace as a matrix.
    """
    E = remove_zero_rows(rref(A, precision=precision), precision=precision)

    n_rows, n_cols = E.shape
    if n_rows == n_cols:
        # E should be the identity, so the nullspace is trivial.
        return []


    # Negative of identity.  We'll replace the rows with rows from E when appropriate.
    I = -np.identity(n_cols)

    # keep track of rows of E
    E_row = 0

    # columns with leading 1's
    lead_cols = [x[1] for x in find_leaing_ones(E)]
    for i, row in enumerate(I):
        if i in lead_cols:
            # if we have leading 1, use E's row
            I[i] = E[E_row]
            E_row += 1

    # set a matrix of right dimensions for the result
    res = np.zeros((n_cols - n_rows, n_cols))

    # select the relevant columns (with "leading -1's")
    for i, col in enumerate([j for j in range(n_cols) if j not in lead_cols]):
        res[i] = I[:, col]

    # make the base "nice"
    return rref(res, precision=precision)

A = np.array([2, 8, 4, 2, 5, 1, 4, 10, -1], dtype=np.float64).reshape(3, 3)
u = np.array([2, 5, 1], dtype=np.float64)

def solve(A, v, precision=1e-9):
    """
    Given a matrix A and vector v (or compatible sizes), solve for A * x = b.
    It gives back a tuple (x, M), where x is a particular solution, and M is a
    matrix that has as rows a basis for the nullspace of A.
    """
    n_rows, n_cols = A.shape

    # create the augmented matrix
    B = np.zeros((n_rows, n_cols + 1))
    B[:, :-1] = A
    B[:, -1] = v.T

    # put it in RREF and  split it back
    E = rref(B, precision=precision)
    AA = E[:, :-1]
    vv = E[:, -1]

    # check for inconsistency
    for row, value in zip(AA, vv):
        # True if zero, False if not
        TF = np.abs(row) < precision
        if np.all(TF) and (np.abs(value) > precision):
            return (None, [])

    # positions of leading ones
    particular_sol = np.zeros(n_cols)
    lead_cols = [x[1] for x in find_leaing_ones(AA, precision=precision)]

    # find particular solution
    for j in range(n_cols):
        if j in lead_cols:
            particular_sol[j] = vv[j]

    # return result (with nullspace)
    return (particular_sol, nullspace(AA, precision=precision))