rref gives an impossible answer

So I am trying to figure out why the following code does give the correct answer. When calculating Q being the matrix P - eye(2), but defined by hand it returns the correct answer however when trying to do the above part, it returns the augmented identity matrix [1 0 0 ; 0 1 0]. I have no idea why it does not return what I expect.

P = [0.9871 0.0027; 0.0129 0.9973];
A = P - eye(2);
B = [A [0;0]];
rref(B)

Q = [-0.0129 0.0027; 0.0129 -0.0027];
rref([Q [0;0]])

Answer 1

[Bruno Luong]

I multiply all your array by 10000, so every inputs are integers

P = [9871 0027; 0129 9973];
A = P - 10000*eye(2);
B = [A [0;0]];
rref(B)

Q = [-0129 0027; 0129 -0027];
rref([Q [0;0]])

ans =

    1.0000   -0.2093         0
         0         0         0


ans =

    1.0000   -0.2093         0
         0         0         0

Now both agree (of course they are now identical when RREF is invoked).

I show you the result are correct by transforming B.

>> C=B/-129

C =

    1.0000   -0.2093         0
   -1.0000    0.2093         0

>> C(2,:)=C(2,:)+C(1,:)

C =

    1.0000   -0.2093         0
         0         0         0

Why you get different results? MATLAB RREF is anything but a serious, it uses rational approximation and fails miserably on floating points. Student is advised to use pencil+paper if not they might fail the exam.

[Sigurd Sundberg]

Does matlab have a tool to perform rref that does not run into problems with floats? Or is it easier to go to integer form to avoid any issues with floats?

The analytical solution(pen + paper) is not an issue at all. Problem here was conflicting outputs, without knowing where they came from.

[Bruno Luong]

What's for? If you want to solve linear system, use "\", if you want to do Reduced row echelon form, use pencil paper.

You can also reduce the toterance

> P = [0.9871 0.0027; 0.0129 0.9973];
A = P - eye(2);
B = [A [0;0]];
rref(B) % Failed

ans =

     1     0     0
     0     1     0

>> rref(B,1e-16) % correct

ans =

   1.000000000000000  -0.209302325581395                   0
                   0                   0                   0

But you would NEVER be sure whereas the result is correct with this toy.

Answer 2

[John D'Errico]

P = [0.9871 0.0027; 0.0129 0.9973];
A = P - eye(2);
A
A =
                   -0.0129                    0.0027
                    0.0129      -0.00270000000000004

Yes. I know that you THINK A is a rank 1 matrix. But is it?

rank(A)
ans =
     2

MATLAB says it is not so.

What did you do? You created a problem where subtractive cancellation amplified irregularities in the least significant bits of your numbers.

Does that happen when you use integers in RREF? No, but only because when you work in integers, you no longer have errors in the least significant bits.

So what does MATLAB do when you type 0.9973?

sprintf('%0.55f',0.9973)
ans =
    '0.9972999999999999642952275280549656599760055541992187500'

It stores the number as a double. Is it EXACTLY 0.9973? NO!!!!!!!! It is the closest approximation to 0.9973 that can be found in a binary storage form. Essentially, in the "number" I write below, if the 1's represent negative powers of 2, the number as it is stored is:

0.11111111010011110000110110000100010011010000000100111

But that is NOT equivalent to 0.9973. It is close. But 0.9973 has no exact finite representation in terms of a series of negative powers of 2, just as 2/3 has no exact finite representation in decimal form.

The other numbers you wrote also have the same failing. For example

sprintf('%0.55f',0.0027)
ans =
    '0.0027000000000000001429412144204889045795425772666931152'

Then you tried to create A, thinking it should be essentially a rank 1 matrix. But it is NOT so. In fact, A is technically invertable. You will get garbage of course.

inv(A)
ans =
     -5.18614284513551e+15     -5.18614284513544e+15
     -2.47782380378693e+16     -2.47782380378694e+16

Essentially it produced the same garbage that any such tool will produce, because while A is not singular, it is also very close to being singular, with errors down in the tiny bits of the numbers, thus: -0.00270000000000004.

What you need to learn is how to work with floating point numbers, but also to understand numerical methods, at least at a basic level. A simple course in numerical analysis would be a great start, then maybe a course in numerical linear algebra. Yes, I know, not gonna happen.

The advice Bruno provided fixes part of this, by removing the floating point errors and the resulting subtractive cancellation errors.

P = [9871 0027; 0129 9973];
A = P - 10000*eye(2);

A
A =
  -129    27
   129   -27

Now A is known EXACTLY, without starting off with crap down in the least significant bits.

Languages like MATLAB can work with integers exactly, because a double exactly represents integers as long as they do not exceed 2^53-1. That is just a fundamental aspect of a double, as defined in the IEEE specs, and many languages use that same storage form, thus the same arithmetic.

rref(A)
ans =
                         1        -0.209302325581395
                         0                         0

So this worked a little better. Not perfectly so.

R = rref(A)
R =
                         1        -0.209302325581395
                         0                         0
                         
Rsym = rref(sym(A))
Rsym =
[ 1, -9/43]
[ 0,     0]

So the EXACT answer is what Rsym provides.

sprintf('%0.55f',R(1,2))
ans =
'-0.2093023255813953598103438480393378995358943939208984375'
    
vpa(sym('-9/43'),55)
ans =
 -0.2093023255813953488372093023255813953488372093023255814

As you can see, they are still just a little different, deviating down around the 17th decimal digit.