Puissance matricielle dans R

Question

En essayant de calculer la puissance d'une matrice dans R, j'ai trouvé que le paquet expm implémente l'opérateur %^%.

Donc x% ^% k calcule la puissance k d'une matrice.

> A<-matrix(c(1,3,0,2,8,4,1,1,1),nrow=3) 

> A %^% 5 
     [,1] [,2] [,3] 
[1,] 6469 18038 2929 
[2,] 21837 60902 9889 
[3,] 10440 29116 4729 

mais, à ma grande surprise:

> A 
    [,1] [,2] [,3] 
[1,] 691 1926 312 
[2,] 2331 6502 1056 
[3,] 1116 3108 505 

en quelque sorte la matrice initiale A a changé à A% ^% 4 !!!

Comment exécutez-vous l'opération d'alimentation matricielle?

Answer 1

J'ai corrigé ce bug dans les sources R-forge (du paquet "expm"), svn rev. 53. ->expm R-forge page Pour une raison quelconque la page Web montre encore rev.52, de sorte que le suivant ne peut pas encore résoudre votre problème (mais devrait dans les 24 heures):

install.packages("expm", repos="http://R-Forge.R-project.org") 

Sinon, obtenir la version svn directement, et installer vous-même:

svn checkout svn://svn.r-forge.r-project.org/svnroot/expm 

Merci à "gd047" qui m'a alerté au problème par e-mail. Notez que R-forge possède également ses propres fonctions de suivi des bogues.
Martint

Answer 2

A^5 = (A^4) * A

Je suppose que la bibliothèque mute la variable d'origine, A, de sorte que la chaque étape consiste à multiplier le résultat-up-à-puis avec la matrice d'origine, R. Le résultat que vous obtenez semble correct, il suffit de les affecter à une nouvelle variable.

Answer 3

Bien que le code source ne soit pas visible dans le paquet, car il est emballé dans un .dll file, je crois que l'algorithme utilisé par le paquet est le fast exponentiation algorithm, que vous pouvez étudier en regardant la fonction appelée matpowfast à la place.

Vous avez besoin de deux variables:

1. result, afin de stocker la sortie,
2. mat, comme une variable intermédiaire.

Pour calculer A^6, depuis 6 = 110 (écriture binaire), à ​​la fin, et result = A^6mat = A^4. C'est la même chose pour A^5.

Vous pouvez facilement vérifier si mat = A^8 lorsque vous essayez de calculer A^n pour 8<n<16. Si oui, vous avez votre explication.

La fonction de package utilise la variable initiale A comme variable intermédiaire mat.

Answer 4

Ce n'est pas une bonne réponse, mais peut-être un bon endroit pour avoir cette discussion et comprendre le fonctionnement interne de R. Ce genre de bug s'est glissé avant dans un autre paquet que j'utilisais.

D'abord, notez que qu'assigner la matrice à une nouvelle première variable ne contribue pas:

> A <- B <-matrix(c(1,3,0,2,8,4,1,1,1),nrow=3) 
> r1 <- A %^% 5 
> A 
    [,1] [,2] [,3] 
[1,] 691 1926 312 
[2,] 2331 6502 1056 
[3,] 1116 3108 505 
> B 
    [,1] [,2] [,3] 
[1,] 691 1926 312 
[2,] 2331 6502 1056 
[3,] 1116 3108 505 

Je suppose que R essaie d'être intelligent qui passe par référence au lieu des valeurs. Pour que cela fonctionne, vous devez faire quelque chose pour différencier A de B:

`%m%` <- function(x, k) { 
    tmp <- x*1 
    res <- tmp%^%k 
    res 
} 
> B <-matrix(c(1,3,0,2,8,4,1,1,1),nrow=3) 
> r2 <- B %m% 5 
> B 
    [,1] [,2] [,3] 
[1,] 1 2 1 
[2,] 3 8 1 
[3,] 0 4 1 

Quelle est la procédure explicite?

Enfin, dans le code C pour le paquet, il y a ce commentaire:

• NB: x sera modifié!L'appelant doit faire une copie si nécessaire

Mais je ne comprends pas pourquoi R permet au code C/Fortran d'avoir des effets secondaires dans l'environnement global.

Answer 5

solution très rapide sans utiliser paquet utilise récursivité: si votre matrice est un

powA = function(n) 
{ 
    if (n==1) return (a) 
    if (n==2) return (a%*%a) 
    if (n>2) return (a%*%powA(n-1)) 
} 

HTH