## Reduced row echelon form technique

Asked by Travis Kocian

### Travis Kocian (view profile)

on 7 Aug 2013
Latest activity Commented on by Lingling Fan

### Lingling Fan (view profile)

on 13 Jun 2019
Accepted Answer by Richard Brown

### Richard Brown (view profile)

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

### Tags

Answer by Richard Brown

### Richard Brown (view profile)

on 7 Aug 2013

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

on 13 Jun 2019
Hi,