Oct 13, 2021

Reduced Row Echelon Form

The blog explains how to find reduced row-echelon form of a matrix in python

Introduction

A matrix being in a row-echelon form means that Gaussian elimination has operated on the rows. A matrix is said to be in column-echelon form if its transpose is in row-echelon form. The specific conditions of RE are:

All rows consisting of only zeros are at the bottom
The leading coefficient (pivot) of a non zero row is always to the right of the leading coefficient of the row above it

RREF (Reduced Row-Echelon Form) has 2 more conditions:

Matrix must be in RE form
All pivots must be 1
Each column containing a leading 1 has zeros everywhere else.

One might achieve the RREF via elementary row operations, and here is a code for that!

Code

Essentially, we want to do this:

def rref(matrix):
  dimension = matrix.shape
  numrows = dimension[0]
  numcols = dimension[1]
  lead = 0 # or pivot
  for r in range(numrows):
      # <END: if you reach last column>

      # <LOOP: to find next pivot>

      # <CONVERT: convert all other rows based on pivot>

Now, lets write code for each one by one. The first part is END, which is essentially that if you dont have any more columns, it means that the matrix is already in a RREF state, hence return.

# END:
if lead >= numcols:
  print(self.matrix)
  return

The loop parts requires you to go through each column in the row to determine if the number is suitable to be a pivot. If you find none, then go to the next column! (this means current column was all zeros and must be put at last). If there are no pivots at all to be found, that means you have already reached RREF

# LOOP:
index = r
while matrix[index][lead] = 0:
  index += 1
  if numrows == index:
      # No pivots in my row, all zeros
      lead += 1
      if lead == numcols:
          # No pivots at all
          print(matrix)
          return

# Swap if I had to go to the next column
matrix[[index, r], :] = matrix[[r, index], :]

I have my pivot now matrix[r][lead], which I need to convert to 1 via elementary row operations. Once that is done, everything else in the column needs to be converted to 0s as well

# CONVERT:

# R = R/k
matrix[r] = matrix[r] / matrix[r][lead]

for i in range(numrows):
  if i != r:
      # I = I - cR
      matrix[i] = matrix[i] - matrix[i][lead]*matrix[r]

🚶Walkthrough:

matrix: array([[ -12,  12,   2],
               [ 12, -12,  -2],
               [  2,  -2, -17]])
------------------------------------------------------------------------

Iteration 1:
lead = 0, i = 0, r = 0

After Swapping:
matrix: array([[ -12,  12,   2],
               [ 12, -12,  -2],
               [ 2,  -2, -17]])

Dividing by -12:
matrix: array([[ 1,  -1,   0],
               [ 12, -12,  -2],
               [ 2,  -2, -17]])

matrix[i] = [ 12 -12  -2] - 12*[ 1 -1  0]
matrix[i] = [ 12 -12  -2] - 12*[ 1 -1  0]

After Iteration 1:
matrix: array([[ 1,  -1,   0],
               [ 0,   0,  -2],
               [ 0,   0, -17]])


------------------------------------------------------------------------

Iteration 2:
lead = 1, i = 1, r = 1

After Swapping:
array([[ 1,  -1,   0],
       [ 0,   0,  -2],
       [ 0,   0, -17]])

Dividing by -2:
matrix: array([[ 1,  -1,   0],
               [ 0,   0,   1],
               [ 0,   0, -17]])

matrix[i] = [ 1 -1  0] - 0*[0 0 1]
matrix[i] = [  0   0 -17] - -17*[0 0 1]

After Second Iteration:
matrix: array([[ 1, -1,  0],
               [ 0,  0,  1],
               [ 0,  0,  0]])
------------------------------------------------------------------------

Third Iteration:

since lead >= numcols, exit

Final Code:

import numpy as np

class RREF:
  def __init__(self, matrix):
      self.matrix = matrix

  def rref(self):
      dimension = self.matrix.shape
      numrows = dimension[0]
      numcols = dimension[1]
      lead = 0
      for r in range(0, numrows):
          if numcols <= lead:
              print(self.matrix)
              return

          index = r
          while self.matrix[index][lead] == 0:
              index += 1
              if numrows == index:
                  index = r
                  lead += 1
                  if lead == numcols:
                      print(self.matrix)
                      return

          # swap if my row is all 0s
          self.matrix[[index, r], :] = self.matrix[[r, index], :]

          # make my first number = 1
          self.matrix[r] = self.matrix[r] / self.matrix[r][lead]

          for i in range(numrows):
              if i != r:
                  print(f"Matrix: self.matrix[i] = {self.matrix[i]} - {self.matrix[i][lead]}*{self.matrix[r]}")
                  self.matrix[i] = self.matrix[i] - self.matrix[i][lead]*self.matrix[r]
          lead += 1

      ic(self.matrix)

m = np.array([
  [-12, 12, 2],
  [12, -12, -2],
  [2, -2, -17]
])
r = RREF(m)
r.rref()

Introduction

Code

🚶Walkthrough:

Final Code:

References: