R: самый быстрый способ настроить матрицу интегралов?

Question

У меня есть функция параметров дерева f(x, y, z) и два ограничения L, U.

Учитывая вектор v, я хочу установить матрицу с элементом M[i, j] = INTEGRAL( f(x, v[i], v[j]) ), где интегралы пределы go от x = L до x = U.

Таким образом, задача состоит из двух элементов:

1. Нам необходимо вычислить интегралы. Мне все равно, как это делается, если это БЫСТРО и достаточно точно. Быстро, быстро, быстро! Какой самый быстрый способ?

2. Нам нужно настроить матрицу M[i, j]. Какой самый быстрый способ?

Пожалуйста, не делайте это вопросом "ВАМ ВАМ ГАССИОНСКУЮ КВАДРАТУРУ или ПРАВИЛА СИМПОНОВ?". Меня не волнует . Скорость здесь важна только. Как бы то ни было, я возьму его, если интегралы по крайней мере с точностью до 1-2 цифр или около того.

Answer 1

Потенциально самое быстрое решение представлено ниже достаточно быстро

Для построения матрицы M я не использовал вложенные for l oop. Идея объясняется в нижней строке, которая, я считаю, значительно ускоряет вычисления

---

Контрольный показатель

Я записал некоторые из возможных решений, и вы можете сравнить их производительность (мое «самое быстрое» решение в method1()).

set.seed(1)
library(pracma)

# dummy data: function f and vector v
f <- function(x) x**3 + cos(x**2)
v <- rnorm(500)

# my "fastest" solution
method1 <- function() {
  m1 <- matrix(0,nrow = length(v),ncol = length(v))
  p <- sapply(seq(length(v)-1), function(k) integral(f,v[k],v[k+1]))
  u <- unlist(sapply(rev(seq_along(p)), function(k) cumsum(tail(p,k))))
  m1[lower.tri(m1)] <- u
  t(m1-t(m1))
}

# faster than brute-force solution
method2 <- function() {
  m2 <- matrix(0,nrow = length(v),ncol = length(v))
  for (i in 1:(length(v)-1)) {
    for (j in i:length(v)) {
      m2[i,j] <- integral(f,v[i],v[j])
    }
  }
  m2 + t(m2)
}

# slowest, brute-force solution
method3 <- function() {
  m3 <- matrix(0,nrow = length(v),ncol = length(v))
  for (i in 1:length(v)) {
    for (j in 1:length(v)) {
      m3[i,j] <- integral(f,v[i],v[j])
    }
  }
  m3
}

# timing for compare
system.time(method1())
system.time(method2())
system.time(method3())

такое, что

> system.time(method1())
   user  system elapsed 
   0.17    0.01    0.19 

> system.time(method2())
   user  system elapsed 
  25.72    0.07   25.81 

> system.time(method3())
   user  system elapsed 
  41.84    0.03   41.89 

---

Принцип

Идея в method1() что вам нужно только вычислить интегралы по интервалам, состоящим из соседних точек в v. Обратите внимание, что интегральные свойства:

• integral(f,v[i],v[j]) равно sum(integral(f,v[i],v[i+1]) + integral(f,v[i+1],v[i+1]) + ... + integral(f,v[j-1],v[j]))

• integral(f,v[j],v[i]) равно -integral(f,v[i],v[j])

В этом смысле, учитывая n <- length(v), вам нужно только выполнять целочисленные операции (что довольно затратно для вычислений по сравнению с транспонированием матрицы или векторным кумулятивным суммированием) n-1 раз (намного меньше, чем choose(n,2) раз в method2() или n**2 раз в method3(), особенно когда n большое).