Matrix power in R

Question

Cercando di calcolare la potenza di una matrice in R, ho trovato che il pacchetto expm implementa l'operatore %^%.

Quindi x% ^% k calcola la potenza k-esima di una 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 

ma, con mia grande sorpresa:

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

in qualche modo la matrice iniziale A è cambiato alla A% ^% 4 !!!

Come si esegue l'operazione di alimentazione della matrice?

Answer 1

Ho corretto quell'errore nei sorgenti R-forge (del pacchetto "expm"), svn rev. 53. ->expm R-forge page Per qualche ragione la pagina web mostra ancora rev.52, quindi quanto segue potrebbe non essere ancora risolvere il tuo problema (ma dovrebbe entro 24 ore):

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

In caso contrario, ottenere la versione svn direttamente, e installare da soli:

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

Grazie a "gd047" chi mi ha avvisato del problema via e-mail. Si noti che R-forge ha anche le proprie funzionalità di tracciamento dei bug.
Martint

Answer 2

A^5 = (A^4) * Un

immagino biblioteca muta variabile originale, A, in modo che il ogni passo consiste nel moltiplicare il risultato-up-till-poi con la matrice originale, R. Il risultato che ritorna sembra corretto, basta assegnarli a una nuova variabile.

Answer 3

Anche se il codice sorgente non è visibile nel pacchetto dal momento che è imballato in una .dll file, credo che l'algoritmo utilizzato dal pacchetto è il fast exponentiation algorithm, che si può studiare, cercando in funzione chiamata matpowfast invece.

avete bisogno di due variabili:

1. result, al fine di memorizzare l'output,
2. mat, come una variabile intermedia.

Per calcolare A^6, poiché 6 = 110 (scrittura binaria), alla fine, result = A^6 e mat = A^4. Questo è lo stesso per A^5.

Si potrebbe facilmente verificare se mat = A^8 quando si tenta di calcolare A^n per qualsiasi 8<n<16. Se è così, hai la tua spiegazione.

La funzione pacchetto utilizza la variabile iniziale A come variabile intermedia mat.

Answer 4

Questa non è una risposta adeguata, ma può essere un buon posto per avere questa discussione e capire il funzionamento interno di R. Questo tipo di bug è stato insinuato prima in un altro pacchetto che stavo usando.

Innanzitutto, nota che semplicemente assegnando la matrice di una nuova variabile prima non aiuta:

> 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 

mia ipotesi è che R sta cercando di essere intelligente passando per riferimento invece di valori. Per farlo funzionare, devi fare qualcosa per differenziare A da 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 

Qual è il modo esplicito per farlo?

Infine, nel codice C per il pacchetto, c'è questo commento:

• NB: x sarà alterato!Il chiamante deve fare una copia se necessario

Ma non capisco perché R lasci il codice C/Fortran avere effetti collaterali nell'ambiente globale.

Answer 5

soluzione molto pratica senza usare alcun pacchetto utilizza ricorsività: se la matrice è una

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

HTH