A^k para la multiplicación de matrices en R?

Question

Supongamos A es una matriz cuadrada. ¿Cómo puedo exponer fácilmente esta matriz en R?

probé dos maneras ya: Ensayo 1 con un corte para-loop y el Ensayo 2 un poco más elegante, pero es todavía muy lejos de Un^k simplicidad.

Ensayo 1

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

Trial 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) 

Answer 1

Aunque Reduce es más elegante, un ciclo for solución es más rápido y parece ser tan rápido como 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 

Por otro lado la solución matrix.power es mucho, mucho más lento:

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

@BenBolker es correcto (una vez más). El for-loop aparece lineal en el tiempo a medida que el exponente aumenta, mientras que la función expm ::% ^% parece ser incluso mejor que el log (exponente).

> 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 

Answer 2

Si es A diagonizable, podría utilizar la descomposición de valores propios:

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)) 
} 

Cuando A no es dia gonalizable, la matriz M (matriz de vectores propios) es singular. Por lo tanto, usarlo con A = matrix(c(0,1,0,0),2,2) daría Error in solve.default(M) : system is computationally singular.

Answer 3

El paquete expm tiene un operador %^%:

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 

Answer 4

Una solución más corto con descomposición de valores propios:

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

Answer 5

De hecho el paquete de la EXPM qué utiliza exponentiation by squaring.

en I pura, esto se puede hacer con bastante eficiencia como tal,

"%^%" <- 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 
} 

Timing esto,

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 

Así, como se esperaba, que son la misma velocidad, ya que utilizan el mismo algoritmo . Parece que la sobrecarga del código r de bucle no marca una diferencia significativa.

Por lo tanto, si no quiere usar expm, y necesita ese rendimiento, entonces puede usar esto, si no le importa mirar el código imperativo.

Answer 6

Aquí hay una solución que no es terriblemente eficiente pero que funciona y es fácil de entender/implementar, funciona en la base y funciona casi al instante a pesar de ser mucho menos eficiente que, p., La solución expm:

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

por ejemplo,

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 

De nuevo, si la eficiencia es de suma importancia, las otras respuestas son muy superiores. Pero esto requiere menos gastos indirectos de codificación y no paquetes para situaciones en las que solo desee calcular algunos poderes de matriz en su trabajo de revisión.