A ^ k pour la multiplication matricielle dans R?

Supposons que Un est certains matrice carrée. Comment puis-je facilement exponentiate cette matrice R?

J'ai essayé de deux façons déjà: Essai 1 avec une boucle for hack et 2 Essai un peu plus élégante, mais elle est encore loin de Unk la simplicité.

Essai 1

set.seed(10)
t(matrix(rnorm(16),ncol=4,nrow=4)) -> a 
for(i in 1:2){a <- a %*% a}

Essai 2

a <- t(matrix(c(0,1,0,0,0,0,1,0,0,0,0,1,0,0,0,0),nrow=4))
i <- diag(4) 
(function(n) {if (n<=1) a else (i+a) %*% Recall(n-1)})(10)

source d'informationauteur hhh | 2012-02-27

  1. 12

    Bien que Reduce est plus élégant, d'une boucle for solution est plus rapide et semble être aussi rapide que expm::%^%

    m1 <- matrix(1:9, 3)
m2 <- matrix(1:9, 3)
m3 <- matrix(1:9, 3)
system.time(replicate(1000, Reduce("%*%" , list(m1,m1,m1) ) ) )
#   user  system elapsed 
#  0.026   0.000   0.037 
mlist <- list(m1,m2,m3)
m0 <- diag(1, nrow=3,ncol=3)
system.time(replicate(1000, for (i in 1:3 ) {m0 <- m0 %*% m1 } ) )
#   user  system elapsed 
#  0.013   0.000   0.014 

library(expm)  # and I think this may be imported with pkg:Matrix
system.time(replicate(1000, m0%^%3))
# user  system elapsed 
#0.011   0.000   0.017

    Sur l'autre côté de la matrice.solution d'alimentation est beaucoup, beaucoup plus lent:

    system.time(replicate(1000, matrix.power(m1, 4)) )
   user  system elapsed 
  0.677   0.013   1.037

    @BenBolker est correct (encore une fois). La boucle for s'affiche linéaire dans le temps que l'exposant augmente alors que la expm::%^% fonction semble être encore mieux que log(exposant).

    > m0 <- diag(1, nrow=3,ncol=3)
> system.time(replicate(1000, for (i in 1:400 ) {m0 <- m0 %*% m1 } ) )
   user  system elapsed 
  0.678   0.037   0.708 
> system.time(replicate(1000, m0%^%400))
   user  system elapsed 
  0.006   0.000   0.006
  2. 12

    La expm paquet a une %^% opérateur:

    library("sos")
findFn("{matrix power}")
install.packages("expm")
library("expm")
?matpow
set.seed(10);t(matrix(rnorm(16),ncol=4,nrow=4))->a
a%^%8
  3. 12

    Si A est diagonizable, vous pouvez utiliser la valeur propre de décomposition:

    matrix.power <- function(A, n) {   # only works for diagonalizable matrices
   e <- eigen(A)
   M <- e$vectors   # matrix for changing basis
   d <- e$values    # eigen values
   return(M %*% diag(d^n) %*% solve(M))
}

    Lorsque A n'est pas diagonalizable, la matrice M (matrice des vecteurs propres) est singulier. Ainsi, en l'utilisant avec A = matrix(c(0,1,0,0),2,2) donnerait Error in solve.default(M) : system is computationally singular.

  4. 2

    Une solution plus courte avec la valeur propre de décomposition:

    "%^%" <- function(S, power)
         with(eigen(S), vectors %*% (values^power * t(vectors)))
  5. 1

    En effet, le expm du package n'utilisez l'exponentiation par la quadrature.

    Dans la pure r, ce qui peut être fait plutôt efficace comme,

    "%^%" <- function(mat,power){
    base = mat
    out = diag(nrow(mat))
    while(power > 1){
        if(power %% 2 == 1){
            out = out %*% base
        }
        base = base %*% base
        power  = power %/% 2
    }
    out %*% base
}

    Délai,

    m0 <- diag(1, nrow=3,ncol=3)
system.time(replicate(10000, m0%^%4000))#expm's %^% function
user  system elapsed 
0.31    0.00    0.31 
system.time(replicate(10000, m0%^%4000))# my %^% function
user  system elapsed 
0.28    0.00    0.28

    Donc, comme prévu, ils sont la même vitesse, car ils utilisent le même algorithme. Il ressemble à la surcharge de la boucle de la r code ne fait pas une différence significative.

    Donc, si vous ne voulez pas utiliser expm, et la nécessité que les performances, alors vous pouvez simplement utiliser, si vous n'avez pas l'esprit en regardant du code impératif.

  6. -1

    Voici une solution qui n'est pas très efficace mais qui fonctionne et qui est facile à comprendre et à mettre en œuvre les œuvres de la base, et fonctionne presque instantanément en dépit d'être beaucoup moins efficace que, par exemple, la expm solution:

    eval(parse(text = paste(rep("A", k), collapse = '%*%')))

    par exemple,

    set.seed(032149)
# force A to be a Markov matrix so it converges
A = matrix(abs(rnorm(100)), 10, 10)
A = A/rowSums(A)
eval(parse(text = paste(rep("A", 1000), collapse = '%*%')))[1:5, 1:5]
#            [,1]       [,2]      [,3]      [,4]       [,5]
# [1,] 0.09768683 0.07976306 0.1074442 0.1284846 0.07789486
# [2,] 0.09768683 0.07976306 0.1074442 0.1284846 0.07789486
# [3,] 0.09768683 0.07976306 0.1074442 0.1284846 0.07789486
# [4,] 0.09768683 0.07976306 0.1074442 0.1284846 0.07789486
# [5,] 0.09768683 0.07976306 0.1074442 0.1284846 0.07789486

    Encore une fois, si l'efficacité est primordiale, les autres réponses sont de loin supérieure. Mais cela nécessite moins de codage frais généraux et pas de paquets pour les situations où vous voulez juste pour calculer la matrice des pouvoirs dans votre brouillon de travail.