更高效地复制一个函数

Question

假设我有以下矩阵。

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函数来处理这些类型的问题。

谢谢。

- shani

---

这是一个巧妙的方法：replicate(5, x <<- x%*%y)。不过我不确定它能快多少。 - James

听起来像是Rcpp的工作。 - baptiste

R 在矩阵运算方面将会比 C++ 更加优越。 - hedgedandlevered

@hedgedandlevered，我喜欢你对R语言的信心。我认为这可能取决于问题的规模、2和5的价值。 - baptiste

Answer 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

- Roland

---

我认为这不正确。如果有什么的话，那就是将 y 进行指数运算。 - James

x 仍然需要先出现（矩阵乘法不满足交换律）。由于这里的 x 和 y 是相同的，因此不会影响答案。 - James

Answer 2

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

将会解决问题。

- Sven Hohenstein

---

问题在于你需要许多 y 的副本，这可能会变得很困难。 - James

@James 同意。这是这种方法的主要缺点。 - Sven Hohenstein

Answer 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

- Roland

---

看起来这只会让事情变得更糟 :). 不过还是给你一个加1，因为你尝试了。 - David Arenburg

@DavidArenburg 我其实很惊讶，使用RcppEigen，我只需要5分钟就能做出一些东西，其速度仅比Martin Maechler合著的代码稍慢。 - Roland

这只能证明你是一位真正的专业程序员。虽然我猜测 %^% 是一个更通用的操作符，它还会进行许多检查以避免潜在的错误等。 - David Arenburg