Параллельное многомерное числовое интегрирование в R?

Question

Я работаю над проблемой, которая требует численного интегрирования двумерной функции, где каждая оценка функции занимает около 1 минуты. Поскольку численное интегрирование на одном ядре будет оценивать функцию от тысячи до десятков тысяч раз, я хотел бы распараллелить вычисления. Сейчас я использую метод bruteforce, который вычисляет наивную сетку точек и складывает их с соответствующими множителями областей. Это определенно неэффективно, и я подозреваю, что любой современный многомерный алгоритм цифровой интеграции сможет достичь такой же точности с гораздо меньшим количеством оценок функций. В R есть много пакетов, которые бы более эффективно и точно вычисляли интеграцию с использованием 2-d (например, R2Cuba), но я не нашел ничего, что можно легко распараллелить в кластере с управляемыми очередями заданий SGE. Поскольку это лишь небольшая часть более крупной проблемы исследования, я хотел бы посмотреть, можно ли это сделать с разумными усилиями, прежде чем я попытаюсь распараллелить один из методов, основанных на кубатурном правиле, в R.

Answer 1

Я обнаружил, что с использованием sparse grid достигается лучший компромисс между скоростью и точностью в многомерной интеграции, и он легко разбивается на кластер, поскольку он не предусматривает каких-либо последовательных шагов. Он не будет таким точным, как другие последовательные алгоритмы интеграции, но он намного лучше, чем наивный метод, поскольку он обеспечивает гораздо более редкую сетку точек для расчета на каждом ядре.

Следующий R-код относится к двумерной интеграции, но может быть легко модифицирован для более высоких измерений. Функция приложения к концу может быть легко распараллелена на кластере.

sg.int<-function(g,...,lower,upper) 

{ require("SparseGrid") 

lower<-floor(lower) 

upper<-ceiling(upper) 

if (any(lower>upper)) stop("lower must be smaller than upper") 

gridss<-as.matrix(expand.grid(seq(lower[1],upper[1]-1,by=1),seq(lower[2],upper[2]-1,by=1))) 

sp.grid <- createIntegrationGrid('KPU', dimension=2, k=5) 

nodes<-gridss[1,]+sp.grid$nodes 

weights<-sp.grid$weights 

for (i in 2:nrow(gridss)) 

{ 
    nodes<-rbind(nodes,gridss[i,]+sp.grid$nodes) 

    weights<-c(weights,sp.grid$weights) 

    } 

    gx.sp <- apply(nodes, 1, g,...) 
    val.sp <- gx.sp %*%weights 
    val.sp 
}