刚才qq群里有人问,R里面有直接算矩阵乘方的函数么?
我居然怎么也没想起来,没辙了只好临时写了一个:
"%^%" <- function(mat, pow) {
if (!is.matrix(mat)) mat <- as.matrix(mat)
stopifnot(!diff(dim(mat)))
if (pow < 0) {
pow <- -pow
mat <- solve(mat)
}
pow <- round(pow)
switch(pow + 1, return(diag(1, nrow(mat))), return(mat))
get.exponents <- function(pow)
if (pow == 0) NULL else c(k <- 2^floor(log2(pow)), get.exponents(pow - k))
ans <- diag(nrow(mat))
dlog2exp <- rev(-diff(c(log2(get.exponents(pow)), 0)))
for (j in 1:length(dlog2exp)) {
if (dlog2exp[j]) for (i in 1:dlog2exp[j]) mat <- mat %*% mat
ans <- ans %*% mat
}
ans
}
小矩阵怎么算都没啥区别,大矩阵可以看出速度差异来。
a = matrix(rnorm(1000*1000),1000)
system.time(a%^%100)
[1] 8.82 0.21 9.24
system.time(a%^%34)
[1] 6.87 0.17 7.19
system.time({
ans <- a
for (i in 1:33) ans <- ans %*% a
})
[1] 31.73 0.65 32.79
以下为同一台机器上MATLAB(默认的优化,我并不保证这种比较是公平的)的计算速度
a = normrnd(0,1,1000,1000);
tic;
b = a^34;
toc
Elapsed time is 13.195601 seconds.
不过我总觉得矩阵乘方这种程序肯定已经有人写过了吧?而且估计应该已经优化得很好了。是不是blas/lapack什么的已经有了?我怎么找不到?