A ^ k для умножения матриц в R?

Question

Предположим, A - это некоторая квадратная матрица.Как я могу легко возвести в степень эту матрицу в R?

Я уже пробовал два способа: Trial 1 с хаком for-loop и Trial 2 немного более элегантно, но это все еще далеко от простоты A ^k .

Пробная версия 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

Если A является диагонализируемым, вы можете использовать разложение по собственным значениям:

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

Когда A не диагонализуем, матрица M (матрица собственных векторов) является сингулярной. Таким образом, использование его с A = matrix(c(0,1,0,0),2,2) даст Error in solve.default(M) : system is computationally singular.

Answer 2

Хотя Reduce более элегантно, решение for-loop быстрее и, кажется, работает быстрее, чем 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 

@ BenBolker верен (еще раз). Цикл for выглядит линейным во времени по мере роста показателя степени, тогда как функция expm ::% ^% оказывается даже лучше, чем log (показатель степени).

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

Действительно, пакет expm использует возведение в степень в квадрате .

В чистом г это можно сделать довольно эффективно, например, так:

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

Итак, как и ожидалось, они имеют одинаковую скорость, потому что используют один и тот же алгоритм. Похоже, что издержки зацикливания кода r не имеют существенного значения.

Так что, если вы не хотите использовать expm и нуждаетесь в такой производительности, тогда вы можете просто использовать это, если не возражаете взглянуть на императивный код.

Answer 5

Более короткое решение с разложением по собственным значениям:

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

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

Опять же, если эффективность имеет первостепенное значение, другие ответы намного выше. Но это требует меньших накладных расходов на кодирование и не требует пакетов для ситуаций, когда вы просто хотите вычислить некоторые возможности матрицы в своей работе.