R-计算外产品矩阵的二次形式的有效方法

Question

我当前正在通过取已知向量并使用两个外部产品的元素乘法作为输入矩阵来计算二次形式。具体来说，我的代码看起来像这样

set.seed(42)  # for sake of reproducibility
library(emulator) 
Fun <- function(a,b) sqrt((1/(2*pi)))*exp(-0.5*(a-b)^2)
n <- 5000
x <- rnorm(n)
y <- rnorm(n)
u <- rnorm(n)
I <- quad.form(outer(x,x,Fun)*outer(y,y,Fun),u)

这很慢，随着n的增加，问题变得明显更糟。据我所知，引起问题的部分是外部(x，x，fun(*外部(y，y，fun(术语。

有什么办法可以加快速度？

Answer 1

您可以通过使用对称性来一半。我发现快速写入RCPP函数最容易：

#include <Rcpp.h>
using namespace Rcpp;
// [[Rcpp::export]]
NumericMatrix foo(NumericVector x, NumericVector y) {
  double n = x.size();
  NumericMatrix M(n, n);
  for(double i = 0; i < n; i++) {
    for(double j = 0; j <= i; j++) {
      M(i,j) = ((1/(2*M_PI))) *
        (exp(-0.5*(pow(x(i)-x(j), 2) + pow(y(i)-y(j), 2))));
      M(j,i) = M(i,j);
    }
  }
  return M;
}

时间：

set.seed(42)  # for sake of reproducibility
library(emulator) 
Fun <- function(a,b) sqrt((1/(2*pi)))*exp(-0.5*(a-b)^2)
n <- 1e4
x <- rnorm(n)
y <- rnorm(n)
u <- rnorm(n)
system.time({
  I <- quad.form(outer(x,x,Fun)*outer(y,y,Fun),u)
})
#       user      system     elapsed 
#      4.287       1.373       5.687 
system.time({
  J <- quad.form(foo(x, y),u)
})
#       user      system     elapsed 
#      2.232       0.168       2.409 
all.equal(I, J)
#[1] TRUE

进一步的改进将需要并行化(或可能使用某些数学(。