Calculating the rank inside a loop

I'm writing a code that calculates Jordan normal form of a matrix.

In this code I take powers of Mk=(M-λI), where M is a given matrix, λ its eigenvalue, I is an identity matrix. Then I define dimensions of kernels of these matrices and storage them. The problem is that MatLab (R2017a) does not calculate the rank of the matrix in each step correctly.

For example, the ranks of two following matrices should be 1 and 2. But MatLab calculates them as 2 and 2.

Mk = -3.0000   -5.0000   -7.0000
    2.0000    4.0000    6.0000
    1.0000    1.0000    1.0000

ans =  2

Mk = -8.0000  -12.0000  -16.0000
    8.0000   12.0000   16.0000
    0.0000    0.0000    0.0000

ans = 2

Here is my code.

M=[0 -5 -7;2 7 6;1 1 4]
n=size(M,1)
e=eig(M)
e=uniquetol(e,0.00001)
for t=1:length(e)
e(t)
c=[1]
draw=[]
for k = 1:n
Mk= (M-e(t)*eye(n))^k
rank(Mk)
draw(k)=n-rank(Mk)
 c=[1 diff(draw)]
if c(end)==0
see there
break
end
end
end

Answer 1

[John D'Errico]

Without even running the code (since I'm not sure what it wants to do), your problem stems from some invalid (although quite common) presumptions when MATLAB displays a matrix as:

Mk = -8.0000 -12.0000 -16.0000

8.0000 12.0000 16.0000

0.0000 0.0000 0.0000

that the numbers are in fact EXACTLY those numbers as shown.

Second, you make the presumption that uniquetol (or any such tool) can produce floating point numbers that will be exactly 0.0001 apart.

Be careful here, because floating point numbers are almost never exactly what they may look to be. NEVER trust a floating point number unless you know enough not to trust what you see. Trust is a difficult thing to come by for floating point numbers.

By the way, a far better way to display floating point numbers is to use the long g option for format.

format long g

e=eig(M)
e =
          4.99999999999999
          2.99999996776984
          3.00000003223016

So M has what appear to be two eigenvalues, 5 and 3, with 3 being an eigenvalue of multiplicity 2. Replicated eigenvalues are often not well estimated however. And the accuracy for a duplicated one tends to be roughly sqrt(eps(e)) where e is the eigenvalue in question. So I would expect the eigenvalue to be not 2.99999996776984 as uniquetol suggests, but a far better estimate might be to form the average of the two miscreant eigenvalues.

mean(e(2:3))
ans =
                         3

As you can see, that would be a much better estimate of the true eigenvalue. Now, let me run your code.

e2 = e(2)
e2 =
          2.99999996776984
k = 1;
Mk = (M-e2*eye(n))^k
Mk =
         -2.99999996776984                        -5                        -7
                         2          4.00000003223016                         6
                         1                         1          1.00000003223016
rank(Mk)
ans =
     2

k = 2;
Mk= (M-e2*eye(n))^k
Mk =
         -8.00000019338094         -12.0000003223016         -16.0000004512222
          8.00000012892063          12.0000002578413          16.0000003867619
      6.44603144195344e-08      6.44603144195344e-08      6.44603155297574e-08
>> rank(Mk)
ans =
     2

So those powers using the inaccurate eigenvalue estimates are increasing in size as k increases. We need not even go very far before things start to fail. Those matrices are NOT exactly what you think they are. There are more digits down in the weeds to worry about.

Rank is telling the truth.

Instead, I'll make one small change, using a more accurate approximation for the miscreant eigenvalue.

e2 = 3
e2 =
     3
k = 1;
Mk= (M-e2*eye(n))^k
Mk =
    -3    -5    -7
     2     4     6
     1     1     1
rank(Mk)
ans =
     2

k = 2;
Mk= (M-e2*eye(n))^k
Mk =
    -8   -12   -16
     8    12    16
     0     0     0
rank(Mk)
ans =
     1

One thing that should be seen as a clue here - when MATLAB displays a number using a decimal point, followed by seros? For example, 12.0000? When that happens, the number is almost never exactly 12.0000, but is rounded purely for display purposes.

Instead, when numbers are displayed as integers, they often are internally represented as integers, or at least they are very close.

mean(e(2:3))
ans =
                         3 % CLOSE? EXACT?

mean(e(2:3)) == 3
ans =
  logical
   0
   
mean(e(2:3)) - 3
ans =
      1.77635683940025e-15

As I said, close, but not exact. Never trust a floating point number really is the best rule, until you get to the point where you know when and where the rules can be broken.

[zdarova_koresh]

Thank you very much for such a long answer. I dealt with the problem by rounding matrices after each step.

[John D'Errico]

Rounding to an integer works. Floating point numbers in double precision can exactly represent integers that do not exceed flintmax, where flintmax is 2^53-1.