Why does a^Inf in matlab returns the 0 matrix?

Question

It seems a^Inf (matrix to the power of infinity) returns a wrong answer:

>> a=[1/2 1/4 1/4; 1/4 1/2 1/4; 1/4 1/4 1/2 ];
>> a^1e99

ans =

    0.3333    0.3333    0.3333
    0.3333    0.3333    0.3333
    0.3333    0.3333    0.3333

>> a^Inf

ans =

     0     0     0
     0     0     0
     0     0     0

What's going on here?

Answer 1

Short answer

Numerical precision issues in the way the matrix power is computed for non-integer exponent.

Long answer

From mpower documentation,

Z = X^y is X to the y power if y is a scalar and X is square. If y is an integer greater than one, the power is computed by repeated squaring. For other values of y the calculation involves eigenvalues and eigenvectors.

a^inf

inf is probably treated as non-integer, and so the method based on eigenvalues and eigenvectors is applied. Namely, the result is computed as

[S,D] = eig(a);
result = S * D^inf * inv(S);

(Quite probably the inverse matrix is not actually computed, but the method is equivalent to this).

For your a, we get

>> a = [1/2 1/4 1/4; 1/4 1/2 1/4; 1/4 1/4 1/2];
>> [S,D] = eig(a)
>> format long
>> D
D =
   0.250000000000000                   0                   0
                   0   0.250000000000000                   0
                   0                   0   1.000000000000000

which looks innocent enough. But wait:

>> D(3,3)-1
ans =
    -3.330669073875470e-16

Since all entries of D have absolute value strictly less than 1, D^inf gives all zeros:

>> D^inf
ans =
     0     0     0
     0     0     0
     0     0     0

and then so does S * D^inf * inv(S), which explains the result for a^inf.

a^1e99

The exponent 1e99 exceeds the maximum integer that can be exactly represented as a double-precision float (which is 2^53), but it is still represented as an integer:

>> mod(1e99,1)
ans =
     0

Thus a^1e99 is computed by the method of repeated squaring. With this method all entries in the result remain close to 0.3333:

>> a^10
ans =
   0.333333969116211   0.333333015441895   0.333333015441895
   0.333333015441895   0.333333969116211   0.333333015441895
   0.333333015441895   0.333333015441895   0.333333969116211
>> a^100
ans =
   0.333333333333333   0.333333333333333   0.333333333333333
   0.333333333333333   0.333333333333333   0.333333333333333
   0.333333333333333   0.333333333333333   0.333333333333333