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
def gauss_method(MATRIX,rescale_leading_entry='Last'):
    """Describe the reduction to echelon form of the given matrix of rationals.

    MATRIX  matrix of rationals   e.g., M = matrix(QQ, [[..], [..], ..])
    rescale_leading_entry='First' make the leading entries to 1's while doing Gaussisan ellimination
    rescale_leading_entry='Last' (Default) make the leading entries to 1's while reducing

    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
        if (rescale_leading_entry == 'First'
           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 
                M.add_multiple_of_row(changed_row,row,factor)
                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
                    M.add_multiple_of_row(row-1-changed_row,row,factor)
                    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]])
gauss_method(M,rescale_leading_entry='First')

And the output looks like this: selection_169

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
def gauss_method(M=random_matrix(QQ,4,algorithm='echelonizable',rank=3),rescale_leading_entry=False):
    """Describe the reduction to echelon form of the given matrix of rationals.

    M  matrix of rationals   e.g., M = matrix(QQ, [[..], [..], ..])
    rescale_leading_entry=False  boolean  make the leading entries to 1's

    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
        if (rescale_leading_entry
           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 
                M.add_multiple_of_row(changed_row,row,factor)
        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: selection_170Or 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 Selection_171.pngAnd 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

Selection_172.png

Do you think this can help students?

Advertisements

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s