Матричная мощность в R

Попытка вычислить мощность матрицы в R, я обнаружил, что пакет expm реализует оператор % ^%.

Итак, x% ^% k вычисляет k-ю степень матрицы.

> 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

но, к моему удивлению,

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

как-то исходная матрица A изменилась на A% ^% 4!!!

Как вы выполняете операцию питания матрицы?

Ответы

Ответ 1

Я исправил эту ошибку в источниках R-forge (пакета expm), svn rev. 53. → expm R-forge page По какой-то причине веб-страница по-прежнему показывает rev.52, поэтому следующее может пока не отображаться решить вашу проблему (но в течение 24 часов):

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

В противном случае получите версию svn напрямую и установите самостоятельно:

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

Благодаря "gd047", который предупредил меня о проблеме по электронной почте. Обратите внимание, что R-forge также имеет свои собственные средства отслеживания ошибок.
Martint

Ответ 2

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

Во-первых, обратите внимание, что просто присваивание матрицы новой переменной сначала не помогает:

> 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

Я предполагаю, что R пытается быть умным, передавая по ссылке вместо значений. Чтобы на самом деле заставить это работать, вам нужно сделать что-то, чтобы отличить A от 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

Каков явный способ сделать это?

Наконец, в коде C для пакета есть этот комментарий:

  • NB: x будет изменен! Вызывающий должен сделать копию при необходимости

Но я не понимаю, почему R позволяет C/Fortran-коду иметь побочные эффекты в глобальной среде.

Ответ 3

Хотя исходный код не отображается в пакете, поскольку он упакован в .dll file, я считаю, что алгоритм, используемый пакетом, это быстрый алгоритм возведения в степень, который вы можете изучить, посмотрев на функцию с именем matpowfast.

Вам нужны две переменные:

  • result, чтобы сохранить вывод,
  • mat, как промежуточная переменная.

Чтобы вычислить A^6, так как 6 = 110 (двоичная запись), в конце концов, result = A^6 и mat = A^4. Это то же самое для A^5.

Вы можете легко проверить, есть ли mat = A^8 при попытке вычислить A^n для любого 8<n<16. Если это так, у вас есть свое объяснение.

Функция пакета использует начальную переменную A как промежуточную переменную mat.

Ответ 4

Очень быстрое решение без использования какого-либо пакета использует рекурсию: если ваша матрица

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

НТН

Ответ 5

A ^ 5 = (A ^ 4) * A

Я предполагаю, что библиотека мутирует исходную переменную A, так что каждый шаг включает в себя умножение результата-up-till-then с исходной матрицей A. Результат, который вы вернетесь, кажется прекрасным, просто назначьте их новому переменная.