MATLAB中矩阵的不精确功率[重复]答案

【问题标题】：Inexact power of matrix in MATLAB [duplicate]MATLAB中矩阵的不精确功率[重复]
【发布时间】：2013-06-23 22:28:03
【问题描述】：

当我很无聊时，我检查了重新分级 MARKOV 链的转移矩阵的平稳定理。所以我定义了一个简单的，例如：

>> T=[0.5 0.5 0; 0.5 0 0.5; 0.2 0.4 0.4];

平稳定理说，如果您将转移矩阵计算为非常高的幂，那么您将得到其主成分位于行的平稳矩阵。那么让我们试试吧：

>> T^1000
ans =

0.4211    0.3158    0.2632
0.4211    0.3158    0.2632
0.4211    0.3158    0.2632

到目前为止一切都很好。继续：

>> T^1000000000
ans =

0.4211    0.3158    0.2632
0.4211    0.3158    0.2632
0.4211    0.3158    0.2632

好的....很好。让我们再多取一个零：

>> T^10000000000

ans =

0.4210    0.3158    0.2632
0.4210    0.3158    0.2632
0.4210    0.3158    0.2632

？？？有些事情发生了变化......让我们尝试更多：

>> T^10000000000000000

ans =

1.0e-03 *

0.5387    0.4040    0.3367
0.5387    0.4040    0.3367
0.5387    0.4040    0.3367

这是怎么回事，连行的总和都不再是1了

 >> T^10000000000000000000

 ans =

 0     0     0
 0     0     0
 0     0     0

啊，它不见了。

我用 R2011a 试过这个。我想在后台有一些奇特的算法，它近似于矩阵的这种高功率。但这怎么会发生呢？哪种算法在此类计算上执行得如此之快，并在这种极端情况下做出这种行为不端？

【问题讨论】：

它预期并且由于二进制系统中的浮点精度算术。换句话说，你不能用有限的位数精确地表示任何十进制数。一个例子：0.3-0.2-0.1 不完全是 0。
将任何涉及浮点数的东西提高到如此高的幂在数值上是疯狂的。你不是在处理符号计算。
考虑到矩阵的大小T、1/10000000000000000（真正出错的幂的倒数）大约是eps，这也许并不奇怪。根据Z = X^y 的帮助（R2012b）“如果y 是大于一的整数，则通过重复平方计算幂。”我认为这仅在大约y = 2^30 的幂内才成立，然后使用另一种方案，可能是特征值分解。
@natan：我看这里没有重复。问题是关于算法的。 @horcher：这很有趣。如果我经常将矩阵取 2 次方，一切似乎都运行良好。事实上，这需要很多时间。因此，这让我更加了解 MATLAB 计算此解决方案的方式。有什么具体的想法吗？
“如果我经常对矩阵进行 2 次幂运算，一切似乎都运行良好。” 你是认真的吗？您的论点是挥手的，请注意2^21 = 2097152，这是浮点精度算术类型问题的常年重复。

标签： matlab matrix markov-chains markov

【解决方案1】：

它可以使用Eigen Decomposition，这样这样的速度是可能的

真正的计算负载是进行这种分解，然后通过计算特征值的幂，可以很容易地计算出幂。这也解释了为什么将计算划分为较小的幂（如 2）需要更多时间。

【讨论】：

【解决方案2】：

一个正如大家所说，原因是浮点精度。解决方案？可变精度算术，在符号数学工具箱中，如果你有的话。只需像这样使用 vpa 初始化矩阵：

T=vpa([0.5 0.5 0; 0.5 0 0.5; 0.2 0.4 0.4],100);

从T^10000000000000000000 给出正确的输出：

[0.42105263157894736842099364708441, 0.31578947368421052631574523531331, 0.26315789473684210526312102942776]
[0.42105263157894736842099364708441, 0.31578947368421052631574523531331, 0.26315789473684210526312102942776]
[0.42105263157894736842099364708441, 0.31578947368421052631574523531331, 0.26315789473684210526312102942776]

【讨论】：