Potencia de matriz en R

Question

Al tratar de calcular la potencia de una matriz en R, encontré que el paquete expm implementa el operador %^%.

Entonces x% ^% k calcula la potencia k-ésima de una matriz.

> 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 

pero, para mi sorpresa:

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

alguna manera la matriz inicial A ha cambiado a un% ^% 4 !!!

¿Cómo se realiza la operación de potencia de la matriz?

Answer 1

He solucionado ese error en las fuentes de R-forge (del paquete "expm"), svn rev. 53. ->expm R-forge page Por alguna razón la página web aún muestra rev.52, por lo que el siguiente todavía no puede resolver su problema (pero debe dentro de las 24 horas):

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

De lo contrario, obtener la versión SVN directamente, e instalar usted mismo:

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

Gracias a "gd047" que me llamó la atención al problema por e-mail. Tenga en cuenta que R-forge también tiene sus propias instalaciones de seguimiento de errores.
Martint

Answer 2

A^5 = (A^4) * A

supongo que la biblioteca muta la variable original, A, de modo que la cada paso implica multiplicando el resultado-up-hasta-entonces con la matriz original, R. El resultado que obtienes parece estar bien, solo asínalos a una nueva variable.

Answer 3

pesar de que el código fuente no es visible en el paquete, ya que se presenta en una .dll file, creo que el algoritmo utilizado por el paquete es el fast exponentiation algorithm, que se puede estudiar examinado la función llamada matpowfast lugar.

Se necesitan dos variables:

1. result, con el fin de almacenar la salida,
2. mat, como una variable intermedia.

Para calcular A^6, ya 6 = 110 (escritura binaria), al final, result = A^6 y mat = A^4. Esto es lo mismo para A^5.

Puede comprobar fácilmente si mat = A^8 cuando intenta calcular A^n para cualquier 8<n<16. Si es así, tienes tu explicación.

La función de paquete utiliza la variable inicial A como la variable intermedia mat.

Answer 4

Esta no es una respuesta adecuada, pero puede ser un buen lugar para tener esta discusión y entender el funcionamiento interno de R. Este tipo de error se ha deslizado antes en otro paquete que estaba usando.

En primer lugar, cabe destacar que simplemente asignar la matriz a una nueva variable primero no ayuda:

> 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 

Mi conjetura es que R está tratando de ser inteligente que pasa por referencia en lugar de valores. Para lograr que esto funcione, debe hacer algo para diferenciar 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 

¿Cuál es la forma explícita de hacerlo?

Por último, en el código C para el paquete, no es este comentario:

• NB: x se verá alterada!La persona que llama debe hacer una copia si es necesario

Pero no entiendo por qué R permite que el código C/Fortran tenga efectos secundarios en el entorno global.

Answer 5

solución muy rápida sin utilizar ningún paquete está utilizando la recursividad: si su matriz es una

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

HTH