Recursive Matrix Row Reduction

Hi,

I'm having trouble debugging a program that I'm writing to reduce a matrix into Reduced Row-Echelon form (RREF). I'm using recursion and Gaussian Elimination to accomplish this task.

The problem is that the last values of the Matrix 'M' in the recursive loop don't seem to carry outside of the recursion; that is, the operations that I perform on them do not seem to stick. This is easier to show with part of the debug output:

    Getting out of the recursion, M is...
        1.0000    0.8333

    Check on m
        1.0000    0.8333

    At the end of the function, M is...
        1.0000    0.8333

    Getting out of the recursion, M is...
         1     3     2
         0     6     5

Full code of the Gaussian Elimination recursive program as below:

    function [M,pivots] = echelon(M)
        %Reduce matrix to reduced row-echelon form
        %ECHELON(M) Returns reduced row-echelon matrix M and pivots vector 'pivots'

        %Check if input is numerical matrix.  If it isn't, give error.
            narginchk(1,1);
            if ~(ismatrix(M) && isnumeric(M))
                error('Not a numeric array');
            end

        %Check if the matrix is empty.  If it is, return
            pivots=[];
            if isempty(M)
                disp('found empty');
                pivots = [];
                return;
            end

        %Find size of Matrix, calculate tolerance level  
            [m,n] = size(M);
            tolerance = length(M)*norm(M)*eps;
        %Find the pivot variable - it is value x, at position p
            [x,p] = max(abs(M(:,1)));

        if x>tolerance
            M = Exchange(M,1,p);

            M = Divide(M,1,M(1,1));

            for i=2:m
                if M(i,1)==0
                    %do nothing
                else
                  Subtract(M,1,i,M(i,1)/M(1,1))
                end
            end

            disp('Recurring now.');
            [~,newpivots] = echelon(M(2:m,2:n));
            disp('Getting out of the recursion, M is...');
            disp(M);

            pivots(1) = 1;
            for z=1:length(newpivots)
                pivots(end+1) = newpivots(z)+1; 
            end

                for z=2:length(pivots)
                    currentpivot = pivots(z);
                    M = Subtract(M,z,1,M(1,currentpivot)/M(z,currentpivot)) 
                end
                  fprintf('Check on m\n');
                    disp(M);

        else
            disp('the else has been called.');
            [~,newpivots] = echelon(M(:,2:n));
            for z=1:length(newpivots)
                pivots(end+1) = newpivots(z) + 1;
            end
        end
        disp('At the end of the function, M is...');
        disp(M);
    end

If you're familiar with Gaussian row reduction, the subtract, exchange, and divide functions used in the program are all elementary row operations (subtracting row i2 by row i1 * a scalar).

Functions as follows:

function returnMatrix = Subtract(M,i1,i2,r)
if ~(nargin == 4)
    error('# of arguments incorrect.  Exiting.');
end
if ~((i1 >= 1) && (i1 <= size(M,1)) && (i2 >= 1) && (i2 <= size(M,1)))
    error('Row i not in bounds of matrix.  Exiting');
end
if r==0
    error('r is zero.  Exiting.');
end
returnMatrix = M;
returnMatrix(i2,:) = returnMatrix(i2,:) - (r*returnMatrix(i1,:));
end

function returnMatrix = Exchange(M,i1,i2)
if ~(nargin == 3)
    error('Incorrect arguments.  Exiting.');
end
returnMatrix = M;
if (i1==round(i1)) && (i2==round(i2)) && (i1 >= 1) && (i2 >= 1) && (i1 <= size(M,1)) && (i2 <= size(M,1))
    returnMatrix([i1 i2],:) = returnMatrix([i2 i1],:);
else
    error('Incorrect arguments entered for row values.  Exiting.');
end
end

function returnMatrix = Divide(M,i,r)
if ~(nargin == 3)
    error('# of arguments incorrect.  Exiting.');
end
if ~((i >= 1) && (i <= size(M,1)))
    error('Row i not in bounds of matrix.  Exiting');
end
if r==0
    error('r value is zero.  Exiting');
end
returnMatrix = M;
returnMatrix(i,:) = returnMatrix(i,:)/r;
end

Sorry, this is a pretty big post, and I've been working on this small error for hours. I was wondering if someone can help me solve the problem I'm having?

Answer 1

[Alex]

Ok, I have your solution.

1. Your subtract function in the beginning did not return the subtraced matrix.

2. You needed to return the previously calculated matrix (as I commented in the first response). However, since the returned matrix is smaller than the current matrix, you have to account for that.

3. The whole pivot point vector is not needed. The current pivot point is calcuated in the beginning of each cycle, meaning that when later on you need to subtract out the other values of the current row, you just use the current row.

Now, in the solution, I assumed that you wanted the generalized format for solving systems of non-zero equations.

I.e. The format of the input M is assumed to be

[matrix A] = [Column solution B] - such that the last column of the input matrix is the vector B.

I'm using an older verion of matlab, so I commented out some of your debugging to make it easier on myself.

if isempty(M)
    disp('found empty');
    %     pivots = [];
    return;
end        %Find size of Matrix, calculate tolerance level
[rows,cols] = size(M);
tolerance = length(M)*norm(M)*eps;
%Find the pivot variable - it is value x, at position p
[x,p] = max(abs(M(:,1)));
if x>tolerance
    M = Exchange(M,1,p);
    M = Divide(M,1,M(1,1));
    for i=2:rows
        if M(i,1)==0
            %do nothing
        else
            M = Subtract(M,1,i,M(i,1)/M(1,1));
        end
    end
    disp('Recurring now.');
    [new_sub_matrix] = echelon(M(2:rows,2:cols));
    disp('Getting out of the recursion, M is...');

      if(~isempty(new_sub_matrix))
          %check to make sure the sub matrix exists - cannot add an empty
          %matrix
          M(2:rows,2:cols) = new_sub_matrix;
      end
      disp(M);

          %we index over colums, but we count based on the # rows
          for z=2:rows
              %         currentpivot = pivots(z);
              subtract_factor = M(1,z);
              M = Subtract(M,z,1,subtract_factor);
          end
  %     end
      fprintf('Check on m\n');
      disp(M);
      % else
      %     disp('the else has been called.');
      %     [new_sub_matrix,newpivots] = echelon(M(:,2:cols));
      %     for z=1:length(newpivots)
      %         pivots(end+1) = newpivots(z) + 1;
      %     end
  end
  disp('At the end of the function, M is...');
  disp(M);

[Christian]

Hi Alex,

I can't thank you enough for taking the time to help me through this problem. I've managed to work through your code to make sure I've completely understood how this works - this assignment has been very beneficial. I modified it a little bit to return the pivots vector, even though it isn't needed, since that was a required output on the assignment sheet.

I still have one small problem. In the solution that you posted, you commented out the last 'else' statement that runs if x < tolerance. With a matrix of all zeros, the program needs to return back the same matrix and an empty pivots vector.

I programmed the last 'else' statement to run when there is no pivot in the first column of M (when x < tolerance). I'd like to use it to solve the dilemma of the zero matrix. Does this make sense?

Thanks again for your time,

Chris

[Christian]

Alex,

It seems that I spoke to soon. I was able to figure out how the pivots vector updated itself in the latter part of the code, and it works fine now! Thanks for all your help.

Chris

[Alex]

A comment to hopefully help you in the future. Don't use the same variable name (like M) throughout a program. If you rename it as you go through, it'll be easier for you to keep track of what should be taking place with it.

Also, what is the pivot vector represent? When I was first looked at your code, the vector was only an index based vector (like [1 2 3 4]). Is the pivot vector actually suppose to be the "solution" vector?

[Christian]

That makes sense, I'll keep that in mind next time.

The pivot vector contains a list of columns that have pivot values in them.

On another note, I seem to be getting a logic error on some test matrices that I've constructed, but not others. Here is an example of a matrice that is producing an incorrect RREF:

x=[0 3 -6 6 4 -5; 3 -7 8 -5 8 9; 3 -9 12 -9 6 15];

In most cases, it is usually only the last column that produces a logic error. Any ideas?

Answer 2

[Alex]

This is a quick look, but I think your command [~,newpivots] = echelon(M(2:m,2:n)); should be [M,newpivots] = echelon(M(2:m,2:n));

M is a matrix, and Matlab does not pass by reference matrix. So, when you are calling your recursive function, a new copy of M is made, within the scope of the new recursive function execution. That new copy must be returned to the calling function if you want the row operations to be kept.

[Christian]

Hi Alex,

I was considering that in my debugging, so I attempted to fix it the way you proposed. Doing so produces an index out of bounds error in the subtraction statements, because it tries to access a matrix M of size 0x1. In addition, Matlab seems to carry all other elements of my matrix M except for the last bit (it seems that the new copy didn't have to returned for the function to work).

However, if this is the only fault that anyone sees, I may have to re-evaluate how I coded the function. I'm still open to more suggestions!

[bym]

maybe looking at this code might help:

http://www.mathworks.com/matlabcentral/answers/13667-solving-system-of-linear-equation-by-gaussian-elimination