Puissance de la matrice en R
en essayant de calculer la puissance d'une matrice dans R, j'ai trouvé que le paquet expm met en œuvre l'opérateur %^% .
donc x % ^ % K calcule la puissance k-th 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
d'une façon ou d'une autre la matrice initiale A a changé à un %^% 4 !!!
comment effectuez-vous l'opération de puissance de la matrice?
5 réponses
j'ai corrigé ce bug dans les sources R-forge (du paquet "expm" ), svn rév. 53. -- > expm R-forge page Pour une raison quelconque, la page Web montre encore Rev. 52, de sorte que les suivants peuvent ne pas encore résoudre votre problème (mais devrait dans les 24 heures):
install.packages("expm", repos="http://R-Forge.R-project.org")
sinon, récupérez la version svn directement, et installez-vous vous-même:
svn checkout svn://svn.r-forge.r-project.org/svnroot/expm
merci à "gd047" qui m'a alerté sur le problème par e-mail.
Notez que r-forge dispose également de ses propres installations de suivi des bogues.
Martint
ce n'est pas une bonne réponse, mais c'est peut-être un bon endroit pour avoir cette discussion et comprendre le fonctionnement interne de R. ce genre de bogue a déjà surgi dans un autre paquet que j'utilisais.
tout d'abord, notez que simplement assigner la matrice à une nouvelle variable d'abord n'aide 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
à mon avis, R essaie d'être intelligent en passant par référence plutôt que par valeurs. En fait pour obtenir que cela fonctionne, vous devez faire quelque chose à 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 façon explicite de faire cela?
enfin, dans le code C pour le paquet, il y a ce commentaire:
- NB: x sera modifié! L'appelant doit faire une copie au besoin
mais je ne comprends pas pourquoi le code C/Fortran a des effets secondaires dans l'environnement global.
bien que le code source ne soit pas visible dans le paquet puisqu'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:
-
result, afin de stocker la sortie, -
mat, en variable intermédiaire.
pour calculer A^6 , depuis 6 = 110 (écriture binaire), à la fin, result = A^6 et mat = 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 n'importe quel 8<n<16 . Si donc, vous avez votre explication.
la fonction package utilise la variable initiale A comme variable intermédiaire. mat .
solution très rapide sans utiliser de colis utilise la 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
A^5 = (A^4)*A
je suppose que la bibliothèque mute la variable d'origine, A, de sorte que chaque étape implique de multiplier le résultat-jusqu'à-puis avec la matrice d'origine, A. le résultat que vous obtenez semblent très bien, il suffit de les affecter à une nouvelle variable.