更有效地复制函数答案

【问题标题】：Replicate a function more efficiently更有效地复制函数
【发布时间】：2014-11-06 14:10:51
【问题描述】：

假设我有以下矩阵。

x <- matrix(seq(1:4), 2, 2)
y <- matrix(seq(1:4), 2, 2)

我想做以下事情。

for(i in 1:5)
{
  x <- x %*% y
}

但是，这是一个简单的示例。我通常有 X 和 Y 的大矩阵，而且 i 也是一个很大的数字。因此，使用 for 循环可能太耗时了。

是否有人知道对这些类型使用 lapply 或 apply 函数。

谢谢。

【问题讨论】：

这是一个 hacky 方式：replicate(5, x <<- x%*%y)。不知道它会快多少。
听起来像是 Rcpp 的工作
R 在矩阵运算方面将优于 C++
@hedgedandlevered 我喜欢你对 R 的信心。我认为这可能取决于问题的大小，2 和 5 的值。

标签： r matrix apply lapply sapply

【解决方案1】：

library(expm)
x %*% (y %^% 5)
#     [,1]  [,2]
#[1,] 5743 12555
#[2,] 8370 18298

基准测试：

set.seed(42)
x <- matrix(rnorm(1e4), 1e2, 1e2)
y <- matrix(rnorm(1e4), 1e2, 1e2)

fun1 <- function(x, y, j) {
  for(i in 1:j)
  {
    x <- x %*% y
  }
  x
}

fun2 <- function(x, y, i) {
  x %*% (y %^% i)
}

fun3 <- function(x, y, i) {
  Reduce("%*%", c(list(x), rep(list(y), i)))
}


library(expm)
all.equal(fun1(x,y,5), fun2(x,y,5))
#[1] TRUE
all.equal(fun1(x,y,5), fun3(x,y,5))
#[1] TRUE

library(microbenchmark)
microbenchmark(fun1(x,y,30), 
               fun2(x,y,30), 
               fun3(x,y,30), times=10)


#Unit: milliseconds
#          expr       min        lq    median        uq        max neval
#fun1(x, y, 30) 21.317917 21.908592 22.103380 22.182989 141.933427    10
#fun2(x, y, 30)  5.899368  6.068441  6.235974  6.345301   6.477417    10
#fun3(x, y, 30) 21.385668 21.896274 22.023001 22.086904  22.269527    10

【讨论】：

我不认为这是正确的。如果有的话，那就是 y 被取幂。
x 仍然需要先来（矩阵乘法不可以交换）。由于 x 和 y 相同，因此不会影响答案。

【解决方案2】：

Reduce("%*%", c(list(x), rep(list(y), 5)))
#      [,1]  [,2]
# [1,] 5743 12555
# [2,] 8370 18298

会成功的。

【讨论】：

麻烦的是你需要很多y的副本，这可能会被证明是禁止的。
@James 同意。这是这种方法的主要缺点。

【解决方案3】：

只是为了好玩，这里有一个使用 RcppEigen 的解决方案：

C++ 代码：

// [[Rcpp::depends(RcppEigen)]]
#include <RcppEigen.h>
using namespace Rcpp;
using   Eigen::Map;
using   Eigen::MatrixXd;
typedef  Map<MatrixXd>  MapMatd;

// [[Rcpp::export]]
NumericMatrix XYpow(NumericMatrix A, NumericMatrix B, const int j) {
   const MapMatd  X(as<MapMatd>(A)), Y(as<MapMatd>(B));
   MatrixXd X1(X);
   for (int i = 0; i < j; ++i) X1 = X1 * Y;
    return wrap(X1);
}

然后在 R 中：

all.equal(fun2(x,y,5), XYpow(x,y,5))
#[1] TRUE

microbenchmark(fun2(x,y,30), 
               XYpow(x,y,30), times=10)
#Unit: milliseconds
#            expr      min       lq   median       uq      max neval
#  fun2(x, y, 30) 5.726292 5.768792 5.948027 6.041340 6.276624    10
# XYpow(x, y, 30) 6.926737 7.032061 7.232238 7.512486 7.617502    10

【讨论】：

这看起来只会让事情变得更糟:)。 (+1) 不错的尝试
@DavidArenburg 让我印象深刻的是，使用 RcppEigen 我可以在 5 分钟内完成一些东西，这比 Martin Maechler 合着的东西要慢一点。
这只能说明你是一个真正的专业程序员。虽然我猜%^% 是一个更通用的运算符，但它也会进行大量检查以避免潜在的错误等。