MATLAB: Reduced row echelon form technique

rref

I am trying to use a code to calculate the reduced row echelon form of a matrix without the function rref. I make a random matrix A and and then make a matrix new_A = (A-lambda*I). The intent is to eventually find the nullspace of new_A without the null function. When comparing the code to the rref command the results match a majority of the time, however on certain occasions there is a discrepancy in the final column. My understanding is that rref must be unique for a given matrix. Any ideas on what could be the source?
clear all
clc
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%





%   MAKES RANDOM MATRIX A
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
m = 3;
casenum = randi(3,1)
if casenum == 1
    A = randn(m);
end
if casenum == 2
    A = i*randn(m);
end
if casenum == 3
    pownum = randi(2,m);
    A = randn(m) + randn(m).*(i.^pownum);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   MAKES MATRIX new_A FROM (A-lambda*I) & CALCULATES MATLAB NULL
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
I = eye(m,m);
mat_eigs = eig(A);
new_A = A - mat_eigs(1)*I  % currently only uses 1st eigenvalue to test 
mat_null = null(new_A);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   FINDS RREF OF new_A
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
mat_rref = rref(new_A)
tol = 1e-12;   
i = 1;
j = 1;
while (i <= m) && (j <= m)
    % Find value and index of largest element in the remainder of column j
    [p,k] = max(abs(new_A(i:m,j)));
    k = k+i-1;
    if (p <= tol)
        % The column is negligible, zero it out
        new_A(i:m,j) = 0;
        j = j + 1;
    else
        % Swap i-th and k-th rows
        new_A([i k],j:m) = new_A([k i],j:m);
        % Divide the pivot row by the pivot element
        Ai = new_A(i,j:m) / new_A(i,j);
        % Subtract multiples of the pivot row from all the other rows
        new_A(:,j:m) = new_A(:,j:m) - new_A(:,j)*Ai;
        new_A(i,j:m) = Ai;
        i = i + 1;
        j = j + 1;
    end
end
code_rref = new_A

Best Answer

It'll be a difference between yours and MATLAB's choice of tolerance. Your matrix should be rank 2, however I notice that MATLAB is often computing the RREF to be the identity.
If you read the docs for rref they specify what they use as tolerance:
max(size(A))*eps *norm(A,inf))
Related Question