# Step-by-step reduction

One of the things that I always tell my students is to check their solutions when they are done solving a problem. That by itself can mean several things, depending on what the problems is. Of course after solving a system of linear equations, one would plug in the solutions into the original system of equations to see if they are satisfied. The harder part is to come up with a way to check if you’ve done the row-reduction correctly or not. One can easily use Sage to see if the reduced row echelon form is computed correctly. If `A` is your matrix, just input `A.rref()` to find the reduced row echelon form of it.

But still sometimes in the case that the answer is wrong, one might get frustrated by not trying to figure out in which step they’ve made an error, or sometimes they might just get stuck. Considering that row-reduction is a time consuming process, it is plausible to assume that some one can easily give up after a few tries. I have found this nice code in Sage that shows step-by-step Gauss-Jordan elimination process, and actually tells you what to do in each step, from here. Then I did a little bit of customization on it. Here is the final version:

```# Naive Gaussian reduction
"""Describe the reduction to echelon form of the given matrix of rationals.

MATRIX  matrix of rationals   e.g., M = matrix(QQ, [[..], [..], ..])

Returns: reduced form.  Side effect: prints steps of reduction.

"""
M = copy(MATRIX)
num_rows=M.nrows()
num_cols=M.ncols()
show(M.apply_map(lambda t:t.full_simplify()))

col = 0   # all cols before this are already done
for row in range(0,num_rows):
# ?Need to swap in a nonzero entry from below
while (col < num_cols
and M[row][col] == 0):
for i in M.nonzero_positions_in_column(col):
if i > row:
print " swap row",row+1,"with row",i+1
M.swap_rows(row,i)
show(M.apply_map(lambda t:t.full_simplify()))
break
else:
col += 1

if col >= num_cols:
break

# Now guaranteed M[row][col] != 0
and M[row][col] != 1):
print " take",1/M[row][col],"times row",row+1
M.rescale_row(row,1/M[row][col])
show(M.apply_map(lambda t:t.full_simplify()))

for changed_row in range(row+1,num_rows):
if M[changed_row][col] != 0:
factor = -1*M[changed_row][col]/M[row][col]
print " take",factor,"times row",row+1,"plus row",changed_row+1
show(M.apply_map(lambda t:t.full_simplify()))
col +=1

print "Above is a row echelon form, let's keep cruising to get the reduced row echelon form:\n"

for i in range(num_rows):
row = num_rows-i-1
if M[row] != 0:
for col in range(num_cols):
if M[row,col] != 0:
if M[row,col] != 1:
print " take",1/M[row][col],"times row",row+1
M.rescale_row(row,1/M[row][col])
show(M.apply_map(lambda t:t.full_simplify()))
break

for changed_row in range(row):
factor = -1 * M[row-1-changed_row,col]
if factor != 0:
print " take", factor,"times row", row+1, "plus row", row-1-changed_row+1
show(M.apply_map(lambda t:t.full_simplify()))
return(M.apply_map(lambda t:t.full_simplify()))
```

And here is a sample run:

```M = matrix(SR,[[3,-1,4,6],[0,1,8,0],[-2,1,0,-4]])
```

And the output looks like this:

Click here to run it in the sage cell server. Note that the code is implemented to work over rational field since it is being used for pedagogical reasons.

In the comments Jason Grout has added an interactive implementation of the first part of the code:

```# Naive Gaussian reduction
@interact
"""Describe the reduction to echelon form of the given matrix of rationals.

M  matrix of rationals   e.g., M = matrix(QQ, [[..], [..], ..])

Returns: None.  Side effect: M is reduced.  Note: this is echelon form,
not reduced echelon form; this routine does not end the same way as does
M.echelon_form().

"""
num_rows=M.nrows()
num_cols=M.ncols()
print M

col = 0   # all cols before this are already done
for row in range(0,num_rows):
# ?Need to swap in a nonzero entry from below
while (col < num_cols                and M[row][col] == 0):              for i in M.nonzero_positions_in_column(col):                 if i > row:
print " swap row",row+1,"with row",i+1
M.swap_rows(row,i)
print M
break
else:
col += 1

if col >= num_cols:
break

# Now guaranteed M[row][col] != 0
and M[row][col] != 1):
print " take",1/M[row][col],"times row",row+1
M.rescale_row(row,1/M[row][col])
print M
change_flag=False
for changed_row in range(row+1,num_rows):
if M[changed_row][col] != 0:
change_flag=True
factor=-1*M[changed_row][col]/M[row][col]
print " take",factor,"times row",row+1,"plus row",changed_row+1
if change_flag:
print M
col +=1
```

One of the uses of this can be to find the inverse of a matrix step by step. For example:

```var('a b c d')
A = matrix([[a,b],[c,d]])
R = gauss_method(A.augment(identity_matrix(A.ncols()),subdivide=True))```

will give you: Or you can run it for a generic matrix of size 3 (for some reason it didn’t factor i, so I used letter k instead:)

```var('a b c d e f g h i j k l m n o p q r s t u v w x y z')
A = matrix([[a,b,c],[d,e,f],[g,h,k]])
R = gauss_method(A.augment(identity_matrix(A.ncols()),subdivide=True))```

The steps are too long, so I’m not going to include a snapshot, if you are interested look at the steps here. But you can check out only the final result by

`view(R.subdivision(0,1))`

which will give you And to check if the denominators and the determinant of the matrix have any relations with each other you can multiply by the determinant and simplify:

`view((det(A)*R.subdivision(0,1)).apply_map(lambda t:t.full_simplify()))`

to get

Do you think this can help students?