Inexact puissance de la matrice dans MATLAB

Question

Comme je m'ennuyais, j'ai vérifié le théorème stationnaire regrading la matrice de transition d'une chaîne MARKOV. Donc, je définissais un simple, .: par exemple

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

Le théorème stationnaire dit, si vous calculez le transitionmatrix à une puissance très élevée, alors vous aurez à la matrice fixe avec ses principaux composants aux lignes. Essayons donc:

>> T^1000 
ans = 

0.4211 0.3158 0.2632 
0.4211 0.3158 0.2632 
0.4211 0.3158 0.2632 

Tout va bien jusqu'ici. Continuons:

>> T^1000000000 
ans = 

0.4211 0.3158 0.2632 
0.4211 0.3158 0.2632 
0.4211 0.3158 0.2632 

OK .... très bien. Prenons juste un zéro de plus:

>> T^10000000000 

ans = 

0.4210 0.3158 0.2632 
0.4210 0.3158 0.2632 
0.4210 0.3158 0.2632 

??? Quelque chose a changé .... Essayons plus:

>> T^10000000000000000 

ans = 

1.0e-03 * 

0.5387 0.4040 0.3367 
0.5387 0.4040 0.3367 
0.5387 0.4040 0.3367 

Ce qui se passe ici, même la somme des lignes n'est plus 1

>> T^10000000000000000000 

ans = 

0  0  0 
0  0  0 
0  0  0 

Aaaand il est parti.

J'ai essayé ceci avec R2011a. Je suppose qu'il y a un algorithme de fantaisie dans le backround, qui se rapproche de cette puissance élevée des matrices. Mais comment cela peut-il arriver? Quel algorithme effectue aussi vite sur de tels calculs, et fait que ce genre de manque se comporte dans des situations extrêmes comme celle-ci?

Answer 1

Il peut utiliser Eigen Decomposition de sorte qu'une telle vitesse est possible

La charge réelle de calcul fait cette décomposition, puis en calculant les pouvoirs des valeurs propres de la puissance peut être facilement calculée. Cela explique également pourquoi le partage du calcul en plus petites puissances comme 2 prend beaucoup plus de temps.

Answer 2

Comme tout le monde le dit, la raison en est la précision en virgule flottante. La solution? Arithmétique à précision variable, dans la boîte à outils mathématique symbolique, si vous l'avez. initialisant simplement votre matrice en utilisant AVP comme ceci:

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

Donne une bonne sortie de T^10000000000000000000:

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