问R:建立积分矩阵的最快方法？
EN

Stack Overflow用户

提问于 2020-04-23 10:46:42

回答 1查看 289关注 0票数 2

我有一个树参数函数f(x, y, z)和两个限制L, U.

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

所以这个问题有两个因素：

，

，我们需要能够计算积分。我不在乎这是怎么做的，只要它快速和合理准确。快，快，快！最快的方法是什么？

，

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

请不要把这个问题写成"dO yOu WaNt GauSsIan QuaDraTure oR SimPsoNs ruLe?“。我不关心。速度是这里唯一相关的东西。只要积分至少精确到1-2位数或其他什么，我就接受它。

matrix

integration

TDSQL-CMySQL免费体验

云原生数据库TDSQL-C MySQL 2核4G 免费体验15天

回答 1

Stack Overflow用户

发布于 2020-04-23 12:16:25

可能最快的解决方案如下所示

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))

关于执行部分所要求的两个要素

--我猜来自package pracma的integral足够快的
来构建矩阵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)时间或method3()中的n**2时间，特别是当n很大时)。