6

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?

0x90
  • 39,472
  • 36
  • 165
  • 245

1 Answers1

11

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
Luis Mendo
  • 110,752
  • 13
  • 76
  • 147
  • 4
    Awesome answer. +1. As an interesting story, this was a question on my Linear Systems Theory test in graduate school. The matrix was invertible with all of its coefficients being less than 1. The answer was to notice that the matrix can be decomposed using eigendecomposition, and the diagonal matrix of the eigenvalues approached the zero matrix as performing a matrix power of a diagonal matrix is simply the power of each element. – rayryeng Jan 30 '18 at 23:51