Reduced row echelon form technique

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

Answer 1

[Richard Brown]

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))

[Lingling Fan]

Hi,

Thanks for your answer!

How could we change the tolerance of rref, as I did not find in docs.

Is there something like "options" like fmincon setup for tolerance?

Thanks in advance!