Эффективный способ вычисления квадратичных форм для матриц внешнего произведения - PullRequest
Новая
       16

Эффективный способ вычисления квадратичных форм для матриц внешнего произведения

2 голосов
/ 26 июня 2019

В настоящее время я вычисляю квадратичную форму, беря известный вектор и используя поэлементное умножение двух внешних произведений в качестве входной матрицы. Если быть точным, мой код выглядит примерно так 

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) член внутри квадратичной формы.

Есть ли способ ускорить это?

1 Ответ

3 голосов
/ 26 июня 2019

Вы можете вдвое сократить время, используя симметрию. Я считаю, что проще всего быстро написать функцию 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

Для дальнейших улучшений потребуется распараллеливание (или, возможно, использование некоторых математических вычислений).

...