How to Implement Back Substitution for Solving Linear Systems in Python?

Question

I'm currently working on implementing a function for back substitution to solve linear systems represented by augmented matrices in Python. I've been studying algorithms for solving linear systems, and I understand the concept of back substitution, but I'm struggling with the implementation details. This function with parameter input is matrix row echelon form 3 x 3 with the augmented matrix to stack with equals in linear system.

# GRADED FUNCTION: back_substitution

def back_substitution(M):
    """
    Perform back substitution on an augmented matrix (with unique solution) in reduced row echelon form to find the solution to the linear system.

    Parameters:
    - M (numpy.array): The augmented matrix in row echelon form with unitary pivots (n x n+1).

    Returns:
    numpy.array: The solution vector of the linear system.
    """
    
    # Make a copy of the input matrix to avoid modifying the original
    M = M.copy()

    # Get the number of rows (and columns) in the matrix of coefficients
    num_rows = len(M)

    ### START CODE HERE ####
    
    # Iterate from bottom to top
    for row in range(num_rows): 
        # Get the substitution row. Remember now you must proceed from bottom to top, so you must start at the last row in the matrix.
        # The last row in matrix is in the index num_rows - 1, then you need to subtract the current index.
        substitution_row = M[row-1]

        # Iterate over the rows above the substitution_row
        for j in range(row + 1, num_rows): 

            # Get the row to be reduced. The indexing here is similar as above, with the row variable replaced by the j variable.
            row_to_reduce = M[j-1]
            # Get the index of the first non-zero element in the substitution row. This values does not depend on j!
            index = get_index_first_non_zero_value_from_row(substitution_row,row)

            # Get the value of the element at the found index
            value = row_to_reduce[index]

            # Perform the back substitution step using the formula row_to_reduce = None
            row_to_reduce = row_to_reduce - value*M[row]

            # Replace the updated row in the matrix, be careful with indexing!
            M[row] = row_to_reduce
            

    ### END CODE HERE ####

     # Extract the solution from the last column
    solution = M[:,-1]
    
    return solution

A = np.array([[1,-1,0.5],[0,1,1], [0,0,1]])
B = np.array([[0.5], [-1], [-1]])
back_substitution(row_echelon_form(A,B))

#Output
[-1,0,-1]

Must output is [1,0,-1]

Answer 1

I recommend not A and B in the same expression, putting them together can lead to a more complex reading of the code/more frequent appearance of unexpected bugs.

A somewhat cleaner solution could be that, the x[i] from the end to the beginning are calculated, each x[i] is equivalent to (b[i] - M[i] @ x) / M[i, i].

def back_substitution(M, b):
    n = b.size
    x = np.zeros_like(b)

    for i in range(n-1, -1, -1):
        x[i] = (b[i] - M[i] @ x) / M[i, i]

    return x

M = np.array([[1, -1, 0.5],[0, 1, 1], [0, 0, 1]])
b = np.array([[0.5], [-1], [-1]])
back_substitution(M,b)

At the end I get the output:

array([[ 1.],
       [ 0.],
       [-1.]])

Answer 2

You gotta watch it carefully. The reason why there is an error is because

"index = get_index_first_non_zero_value_from_row(substitution_row,row)".

the first parameter of get_index_first_non_zero_value_from_row is n x n matrix not substitution_row.

I know this problem is from Linear Algebra for Machine Learning and Data Science in Coursera. I wish you solve the problem as you want :D