आर में मैट्रिक्स गुणा के लिए ए^के?

Question

मान लीजिए ए कुछ वर्ग मैट्रिक्स है। मैं आर में इस मैट्रिक्स को आसानी से कैसे बढ़ा सकता हूं? एक के लिए लूप हैक और परीक्षण 2 में थोड़ा और अधिक सुंदर ढंग से साथ ट्रायल 1 लेकिन यह अभी भी एक^{कश्मीर} सादगी से एकदम अलग है:

मैं पहले से ही दो तरह से कोशिश की।

ट्रायल 1

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

ट्रायल 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

हालांकि Reduce, और अधिक सुरुचिपूर्ण है के लिए लूप समाधान तेजी से होता है और लगता है एक के रूप में तेजी से होने के लिए के रूप में 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 

दूसरी ओर matrix.power समाधान है बहुत, बहुत धीमी:

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

@ बेनबॉल्कर सही है (फिर भी)। फॉर-लूप समय में रैखिक दिखाई देता है क्योंकि एक्सपोनेंट बढ़ता है जबकि एक्सएमएम ::% ^% फ़ंक्शन लॉग (एक्सपोनेंट) से भी बेहतर होता है।

> 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

तो A है diagonizable, आप eigenvalue अपघटन इस्तेमाल कर सकते हैं:

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

जब एक व्यास नहीं है gonalizable, मैट्रिक्स M (eigenvectors के मैट्रिक्स) एकवचन है। इस प्रकार, A = matrix(c(0,1,0,0),2,2) के साथ इसका उपयोग Error in solve.default(M) : system is computationally singular देगा।

Answer 3

expm पैकेज एक %^% ऑपरेटर है:

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

eigenvalue अपघटन के साथ एक छोटा समाधान:

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

Answer 5

दरअसल expm के पैकेज exponentiation by squaring का उपयोग करता है।

शुद्ध आर में, यह नहीं बल्कि कुशलता से इतने तरह किया जा सकता है,

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

समय इस,

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 

तो, उम्मीद के रूप में, वे एक ही गति कर रहे हैं क्योंकि वे एक ही एल्गोरिथ्म का उपयोग । ऐसा लगता है कि लूपिंग आर कोड के ऊपरी हिस्से में कोई महत्वपूर्ण अंतर नहीं होता है।

तो, यदि आप एक्सएमएम का उपयोग नहीं करना चाहते हैं, और उस प्रदर्शन की आवश्यकता है, तो आप केवल इसका उपयोग कर सकते हैं, अगर आपको अनिवार्य कोड को देखने में कोई फर्क नहीं पड़ता है।

Answer 6

यहां एक ऐसा समाधान है जो बहुत कुशल नहीं है लेकिन काम करता है और समझने/कार्यान्वित करने में आसान है, आधार पर काम करता है, और लगभग कम कुशल होने के बावजूद लगभग तुरंत काम करता है, उदाहरण के लिए, expm समाधान:

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

जैसे

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 

फिर, अगर दक्षता सर्वोपरि है, अन्य उत्तर कहीं बेहतर कर रहे हैं। लेकिन इसके लिए कम कोडिंग ओवरहेड की आवश्यकता होती है और परिस्थितियों के लिए कोई पैकेज नहीं होता है जब आप अपने स्क्रैच काम में कुछ मैट्रिक्स शक्तियों की गणना करना चाहते हैं।