为什么matlab中的^ Inf返回0矩阵？

Question

似乎a^Inf(矩阵与无限的力量)返回错误的答案:

>> 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

这里发生了什么？

Answer 1

简短的回答

计算非整数指数的矩阵幂的方式中的数值精度问题.

答案很长

从mpower文档,

Z = X^y是X的y权力,如果y是标量,X是方形的.如果y 是大于1的整数,则通过重复平方来计算功率.对于其他y计算值,涉及特征值和特征向量.

a^inf

inf可能被视为非整数,因此应用基于特征值和特征向量的方法.即,结果计算为

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

(很可能逆矩阵实际上并未计算,但该方法与此相当).

为了你a,我们得到

>> 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

看起来很无辜.可是等等:

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

由于所有条目的D绝对值都严格小于1,所以D^inf给出全零:

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

然后这样做S * D^inf * inv(S),这解释了结果a^inf.

a^1e99

指数1e99超过可以精确地表示为一个双精度浮点(其中最大整数是 2^53),但它仍然是表示为一个整数:

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

因此a^1e99通过重复平方的方法来计算.使用此方法,结果中的所有条目都保持接近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