R：建立积分矩阵的最快方法？

Posted 2023-03-24

【中文标题】R：建立积分矩阵的最快方法？

【英文标题】：R: fastest way to set up matrix of integrals?

【发布时间】：2020-08-07 05:21:49

【问题描述】：

我有一个树形参数函数f(x, y, z)，还有两个限制L, U。

给定一个向量v，我想用元素M[i, j] = INTEGRAL( f(x, v[i], v[j]) ) 建立一个矩阵，其中积分限制从x = L 到x = U。

所以问题有两个要素：

我们需要能够计算积分。我不在乎这是如何完成的，只要它快速且相当准确即可。快，快，快！！最快的方法是什么？

我们需要设置矩阵M[i, j]。最快的方法是什么？

请不要把这个问题当作“你想要高斯正交还是辛普森法则？”的问题。我不关心。速度是这里唯一相关的东西。不管哪个更快，我都会接受，只要积分至少精确到 1-2 位或其他数字即可。

【问题讨论】：

从表面上看，这可以在双 for 循环中使用 R stats 基本函数 integrate 来完成。只有当我们能看到您正在尝试做的事情的示例时，效率才能真正得到解决。

如果你想要的话，看看我的回答

【参考方案1】：

下面给出了一个可能最快的解决方案

library(pracma)
M <- 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))))
M[lower.tri(M)] <- u
M <- t(M-t(M))

---

关于 OP 要求的两个元素

我猜 integral 来自包 pracma 已经足够快了 为了构建矩阵M，我没有使用嵌套for 循环。这个想法在底线得到了解释，我相信这会显着加快计算速度

---

基准测试

我写下了一些可能的解决方案，您可以比较它们的性能（我的“最快”解决方案在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 次（远低于method2() 中的choose(n,2) 次或n**2 倍于method3()，尤其是当n 很大时。

【讨论】：

我对这个问题的解释有点不同。 f(x,y,z) 具有三个独立的参数。积分超过了固定限制 L 和 U 之间的 x。这种函数的示例是 cos(y*x)*cos(z*x)，积分限制为 0 到 n*pi。 v 是频率向量。我认为您的 method1 不适用于这种解释。

@WaltS 是的，也许你的解释是对的！ OP 没有显示任何问题描述示例。如果情况如您所述，那么嵌套的for 循环可能是不可避免的