CSC 160: Gaussian Reduction Assignment

1. Non-pivoting Gaussian Elimination

This task is basically like Attaway

Write a simple function gauss_reduce() with prototype

function solution_vec = gauss_reduce(param_mat, const_vec)
to solve a system of linear equations using direct Gaussian reduction. This function takes a square matrix of parameters, and a vector of constants, and returns a solution vector. (There are other possibilies, e.g. a function that takes a single matrix representing both the coefficients and the constant vector, but let's not worry about that for now). Your function should be able to take systems of any size (up to practical limits). It should check that the matrix is square, and that the constant vector has the appropriate number of elements, printing an error message and returning the zero vector if these solvability conditions are violated. Your program should also print an error message and return the zero if it runs into a zero pivot which might indicate a singular system, though does not necessarily imply that.

In general zero pivots may not be obvious, given the initial coefficient matrix. If the top left element, the first pivot, is zero that IS obvious, but the process of reducing the first row may set the (2,2) element (the second pivot) to zero when it was non-zero before, and so on. The potential pivot element must be checked as the first step in the column-reduction operation (the reduce_column() function described below).

What to turn in for this part:

Use rand and ones to make a parameter matrix A, a 4x4 random matrix with elements between -1 and 1. It turns out that it's really unlikely these matrices will cause any zero-pivot problems. Also make constant vector B, a 4x1 similarly-random column vector. Run your gaussian-elimination solver on the (A, B) system to get your solution vector X. Check X two ways:

  1. Check that AX = B.
  2. Use the Matlab built-in linear equation solver X = A\B to see if Matlab gets it right.

Why write your own solver if one exists already? First, it's a well-known coding exercise at about the right level (it's assigned in the text, for instance). More importantly, you have actually written a routine that's not built-in, and that's a non-pivoting Gaussian elimination routine. "Normally" we always want to pivot, but the extra credit part of this assignment (below) tries to establish exactly why, since there seem still to be open questions about it. It's really not unusual to have to build your own version of commonly-available functions just so you can vary, instrument (collect statistics), and experiment with them.

  • Again, use rand and ones to make a parameter matrix A, a 4x4 random matrix with elements between -1 and 1. Make constant vector B, a 4x1 similarly-random column vector. Set
    A(1,1) = 0. Run your gaussian-elimination solver on the (A, B) system. You should get an error for zero pivot.

    Now swap the first two rows of A and B to get Aswap and Bswap and try your gaussian-elimination program on the Aswap and Bswap system. It should run OK. Recall that row-swapping is an Elementary Row Operation (ERO) that simply re-arranges the order of the same equations, so the solution to the system is the same as the unswapped system. Congratulations, you've gotten over a zero-pivot problem. Remember the answer (X) vector.

    One bad thing about built-in Matlab operations is we have no idea what's going on inside. One good thing is that they usually are reliable and "have seen it all" and are prepared for nasty inputs. Since X = A-1B, in some metaphorical sense X = B/A and in Matlab, writing the last equation with a "divided into (\)" rather than a "divided by (/)" operator invokes linear equation solution.

    So, try Xswap = Aswap\Bswap;, which should work since you fixed the zero-pivot problem. Then similarly try X = A\B. What can you conclude about the built-in \ operator?

    In fact, your row swap manually did one partial pivot operation on the first column of A. Partial pivoting is the use of row swaps to avoid zero pivots. Before eliminating the below-diagonal elements of any column starts, partial pivoting usually swaps the maximum column element in or below the main diagonal of the column up to the main diagonal to act as the pivot. The next (extra credit) part of the assignment has you implement partial pivoting and investigate whether it improves, as is often claimed but not proved, the numerical accuracy of solutions.

    Coding Hints

    Your gaussian-elimination function should use secondary functions called by your main function as a structuring mechanism. A good place to start is to note that after the initial parameter checking, the process consists of two main steps. First, reduce the coefficient matrix to upper triangular form (modifying the constant vector in parallel). Second, perform the back-substitution to obtain the solution vector from the upper triangular matrix and the (modified) constant vector. We can write a secondary function to perform each step. Our main function thus starts out looking like this:

    function solution_vec = gauss_reduce(param_mat, const_vec)
      % check for consistent size of param_mat and const_vec
      % reduce coefficient matrix to upper triangular form, modifying the
      % constant vector appropriatly
      [ut_mat, new_const_vec] = ut_reduce(param_mat, const_vec);
      % Compute the solution vector using back substitution
      solution_vec = back_subst(ut_mat, new_const_vec);
      % we are done

    The ut_reduce() function uses its own subsidiary functions. Specifically, you should write a function called reduce_column() with prototype

    function [new_mat, new_const_vec] =
        reduce_column(param_mat, const_vec, column)
    that returns a modified matrix (and constant vector) in which the input matrix has been reduced so that all the elements below the (col,col) diagonal element are zero. Initially it checks to see the (col, col) element is NOT zero, of course! This can be used iteratively to modify the original matrix and constant vector to produce an upper triangular form. Note that using a function to modify the matrix and vector passed as parameters in the following way,
    [cur_mat, cur_vec] = reduce_column(cur_mat, cur_vec, cur_col);
    is perfectly legitimate (as long as you don't need the old partial solutions) and a good way of structuring the process.

    The reduce_column() function might itself call yet another function to reduce a specified column in a row to 0 by adding a multiple of another row, returning again, both a modified coefficient matrix and a modified constant vector. The prototype would look like

    function [new_mat, new_const_vec] =
        reduce_row_at_col(param_mat, const_vec, ...
                          col, row_added, row_reduced);
    This function adds a multiple of row_added to row_reduced so that the specified col in row_reduced is 0. The same operation is carried out on the corresponding position of the constant vector. (Basically like the "ERO" Attaway describes at the bottom of page 343)

    The back-substitution function can similarly be constructed using a subsidiary function that takes the upper triangular matrix, the (modified) constant vector, and a partially-filled-in-from-the-end solution vector, (entries from col + 1 to the end already known) and produces a modified partial solution vector with entry col now filled in. The prototype would look like

    function new_part_solution_vec = 
        back_subst_for_col(ut__mat, new_const_vec,...
                            column, part_solution_vec)
    Note that if you call for a column to be filled without previously filling in all the higher columns, the function will probably not work as desired, so you need to be a little careful how you use it. Since it is your function, and not one you are publishing to the world, this is OK. Even so, you should leave yourself a note in comments.
  • Despite all the above verbiage, the amount of programming needed is small. CB did down through Dessert (all of Main Course plus partial pivoting, not the experiments) in 40 lines (32 for elimation with pivoting, 8 for back-substitution). He didn't use exactly these functions, but they are similar in spirit and small. Most are 4-liners; there's a 7 and a couple of 6's, longest is 8 for back-sub. This doesn't count his dense and helpful comment lines of course. There are some terse but not-well-structured solutions on the web too.

    Extra Credit:

    Partial Pivoting

    Write a second version of your function called gauss_reduce_pp() that uses partial pivoting. Write a secondary function that identifies and swaps in the correct pivot row for a specified working column (modifying a matrix and constant vector) to help you do this. The prototype would be something like
    function [new_mat, new_const_vec] = pivot(param_mat, const_vec, col);
    A function to swap a specified pair of rows (modifying both the coefficient matrix and constant vector) with prototype
    function [new_mat, new_vec] = swap_rows(mat, vec, row1, row2);
    will be useful.

    Also, write a function with prototype

    function [param_mat, const_vec] = random_test_case(n)
    to generate random test cases of a specified size n, specifically coefficient matrices and constant vectors with values that are floating point values (not integers) between -100.0 and 100.0.

    Correctness, Timing and Distribution Tests

    Now that you've got your basic tools, the research begins.

    1. Run your programs on examples of 2x2, 3x3, and 4x4 systems, and check that the solutions are correct with a function that takes the coefficient matrix, the solution vector, and the constant vector, multiplies the coefficients by the solution and checks that the result is equal to the constant vector. You may want NOT to use a straight == equality check, but a check using eps() (see help eps.) The issue is that floating point numbers may be really close (effectively equal) but not exactly the same, due to having been calculated different ways with roundoff error affecting the intermediate and final results.
    2. Use tic, toc to time your system running on 10x10, 100x100, and 1000x1000 systems. Do the results seem to support the theoretical result that the running time on an NxN system is proportional to N3? A plot over a range of sizes could be more convincing. Repeat using the partial pivoting program and compare the times. What percentage overhead does adding pivoting entail?
    3. Run your solver on 100 random 3x3 examples (generated by your random_test_case() function), and compute the mean and standard deviation of each of the 3 solution variables. That's three mean values and three stds all told. (Question: What should the means be like (remember the characteristics of the inputs!)). (Hint: They should be similar, and the mean should be small compared to 100.0). Presumably you will use a program (function or script) to do this. Otherwise you will be typing till the cows come home...

      Of course you only care about statistics over those 100-long trials. For each size of perturbation you're presumably interested in the difference pivoting makes to the accuracy of the answers. So we'd expect a total of maybe 20 numbers to be your final output: something like the mean and standard deviation of the absolute value of the 100 errors for each of the 5 perturbation sizes and two pivoting conditions.

      Twenty numbers still sounds like a lot, so a good idea is to make one plot with two data series plotted. One is means and std error-bars (matlab can plot those, use help or the documentation) for pivoting. The other gives the values for the no-pivoting case. So we'll see a single visual with all the information, readily comparable...each "data point" represents 100 experiments or so.

    Perturbation Experiment and Effect of Pivoting

    1. For each of 100 random 10 x 10 systems generated by your program (i.e. with coefficient values between -100.0 and 100.0), do the following. First find the solution to the system and record it. Then perturb the system by adding a (different) random (real-valued) number between -1.0 and 1.0 to each of the matrix coefficients and each constant element. (Note that this range is 1% of the original range). Solve the perturbed system. Compute the distance between the solution to the perturbed problem and the solution to the original problem , i.e. compute
      E = √ (Δ x 1 2 + Δ x 2 2 + Δ x 3 2 + ... + Δ x 10 2 ).
      Here Δ x is the vector difference between the solution vectors of the original and perturbed problems. This square root of the sum of its squared elements (E above) is a measure of the error. Compute the mean and standard deviation of these 100 errors (two numbers: one mean, one standard deviation).

      Now we think it's pretty likely that you won't get a significant result (interesting differences in the mean and stdevs) with that perturbation in the range [-1.0, 1.0], so if not try smaller perturbations -- range [-0.01, 0.01]? [-0.00001, 0.00001]? etc. to see if you can uncover anything with them. If you're getting into the 160 spirit, you'll recognize a superb opportunity to exploit this idea and just go ahead and run a bunch of perturbations, each smaller by a factor of 10, and make a semilog plot of your results. This is not only interesting, but would add substantial visual 'reader appeal' (or 'grader appeal, if you care about that) to your writeup.

      A small aside. Make sure that an error exit from your simple program because it encountered a zero pivot (or other reason) does not contaminate the data. The probability of a zero pivot occuring is quite low, but if it did, the zero return vector would bollix the mean. Re-run the experiment (with new random values) if this ever happens. Actually, if you find this occuring on random inputs, check your programs (and random data generation) for errors, as the probability is really extremely low.

    2. Repeat the previous perturbation experiment using the program with partial pivoting. The random systems do not have to be the same ones used above, although they can be. Is there any significant reduction in the mean error? The idea of this experiment is to test the (oft heard) claim that using partial pivoting reduces the sensitivity of the solution process to error.

    Ill-Conditioned Systems and More

    1. A random ill-conditioned 10 x 10 system can be produced by first generating the a coefficient matrix (and corresponding constant vector) at random, using values between -100.0 and 100.0 as before. Then randomly select 3 of the 10 rows, (they all have to be different) and make one of them almost the average of the other two by taking the sum/2 (average) plus a random value between -1.0 and 1.0 (1% of the original range). How to pick 3 random positions from 10 is a little bit of a puzzle, but not that hard (first pick one of 10, then one of the 9 remaining, then ...) Repeat the perturbation experiment using 100 such ill-conditioned systems. (that is compare answers to original and perturbed systems using the ill conditioned systems you have made, rather than purely random ones)

      Test both the simple program and the version using partial pivoting. How do the mean and standard deviation of the error compare to that of the (mostly) well conditioned random systems? How does the effect of partial pivoting compare? Again, if a random 'noise' addition in the 1% range doesn't yield any discernable differences in behavior, try a smaller range... 0.01% ? 0.00001%? etc. No guarantees, but it IS research and it DOES give you more to write about and it's easy to do. Again again, using a bigger set of random ranges, probably multiples of 10 as above, is a chance to show off your plotting prowess and make a scientific point visually, always better than a table or a bunch of numbers stuck in a paragraph of text.

      For more, repeat the pertubation experiment (the first one, not the one above using ill-conditioned matrices) using 100 randomly generated systems of size 5, 10, 20, 50, 100, [200, 500, 1000] (the last 3 might take too much time, depending on your implementation). Is the mean error correlated with the size of the system? Attempt to explain.

    What to Hand In

    See the Universal Hand-In Guide.

    Briefly, for the Code component of this assignment, you'll need a .zip (not .rar) archive with code files and README.

    Submit the writeup as a single, non-zipped .PDF file.

    Submit before the drop-dead date for any credit and before the due date for partial-to-full (or extra!) credit.

    Check immediately to see BB got what you sent!

    Last Change: 5/6/2011: CB