Мощность матрицы в R

Question

Пытаясь вычислить мощность матрицы в 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 !!!

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

Answer 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

Answer 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 для пакета, есть этот комментарий:

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

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

Answer 3

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

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

НТН

Answer 4

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

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

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

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

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

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

Answer 5

Неэффективная версия (поскольку более эффективно сначала диагонализировать матрицу) в base без особых усилий:

pow = function(x, n) Reduce(`%*%`, replicate(n, x, simplify = FALSE))

Я знаю, что этот вопрос конкретно о старой ошибке в expm,но в настоящее время это один из первых результатов для "матрицы мощности R", поэтому, надеюсь, эта небольшая стенография может быть полезна для кого-то еще, кто в итоге просто ищет быстрый способ запустить матрицу мощности без установки каких-либо пакетов.

Answer 6

A ^ 5 = (A ^ 4) * A

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