Matriz de potencia en R

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.

Primero, tenga en cuenta 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 pasando por referencia en lugar de valores. Para que esto funcione, necesitas 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

¿Qué es el ¿manera explícita de hacer esto?

Finalmente, en el código C para el paquete, hay este comentario:

• NB: x será alterado! 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.

Aunque el código fuente no es visible en el paquete ya que está empaquetado en un .dll file, creo que el algoritmo utilizado por el paquete es el algoritmo de exponenciación rápida, que puede estudiar mirando la función llamada matpowfast en su lugar.

Necesitas dos variables:

1. result, para almacenar la salida,
2. mat, como variable intermedia.

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

Usted podría 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 package utiliza la variable inicial A como variable intermedia mat.

Solución muy rápida sin usar ningún paquete está usando recursividad: si su matriz es 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

Supongo que la biblioteca muta la variable original, A, de modo que cada paso implica multiplicar el resultado-hasta-entonces con la matriz original, A. El resultado que obtiene parece estar bien, solo asígnelos a una nueva variable.