Попытка вычислить мощность матрицы в 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!!!
Как вы выполняете операцию питания матрицы?
Я исправил эту ошибку в источниках 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
Это не правильный ответ, но может быть хорошим местом, чтобы обсудить это и понять внутреннюю работу 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 для пакета есть этот комментарий:
Но я не понимаю, почему R позволяет C/Fortran-коду иметь побочные эффекты в глобальной среде.
Хотя исходный код не отображается в пакете, поскольку он упакован в .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.
Очень быстрое решение без использования какого-либо пакета использует рекурсию: если ваша матрица
powA = function(n)
{
if (n==1) return (a)
if (n==2) return (a%*%a)
if (n>2) return ( a%*%powA(n-1))
}
НТН
A ^ 5 = (A ^ 4) * A
Я предполагаю, что библиотека мутирует исходную переменную A, так что каждый шаг включает в себя умножение результата-up-till-then с исходной матрицей A. Результат, который вы вернетесь, кажется прекрасным, просто назначьте их новому переменная.