Среднее значение столбца data.table, указанное с использованием матрицы
У меня есть data.table, содержащий значения x,y,z в 10000 точек (для этого примера) в единичном кубе, и каждая точка имеет соответствующий атрибут (называемый P). Я использовал nn2 от RANN пакет для поиска индексов k-соседей (до 50) каждой точки в радиусе 0,075 единиц от исходного data.frame (который возвращается в виде матрицы).
library(RANN)
library(data.table)
set.seed(1L) # for reproducible data
DATA <- data.table(runif(10000, 0,1),
runif(10000, 0,1),
runif(10000, 0,1),
runif(10000, 10,30))
colnames(DATA)<-c("x","y","z","P")
nn.idx <- nn2(DATA[,1:3], DATA[,1:3], k=50,
treetype = "kd", searchtype = "radius",
radius = 0.075)$nn.idx
Следующие for Цикл делает работу, но мне было интересно, есть ли способ ускорить это, векторизовав это, поскольку это не будет масштабироваться при применении к> миллионам точек? Проще говоря, я хочу использовать nn.idx получить соответствующий P значения из DATA и рассчитать среднее P который затем присваивается новому столбцу в DATA называется mean.P
for(index in 1:nrow(DATA))
DATA$mean.P[index]<-mean(DATA[nn.idx[index,], P])
В иллюстративных целях следующий код иллюстрирует то, что я пытаюсь вычислить - для всех точек (серые точки) вычислите среднее значение для всех точек (оранжевые + красные точки) в сфере вокруг данной точки (красная точка) и назначьте это к этому моменту (красная точка). Итерируйте по всем точкам, но делайте это эффективным способом, который будет масштабироваться для больших наборов данных.
library(rgl)
rgl.open()
rgl.points(DATA[1500,1], DATA[1500,2], DATA[1500,3], color ="red")
rgl.points(DATA[nn.idx[1500,],1:3], color ="orange", add=T)
rgl.points(DATA[,1:3], color ="lightgray", alpha=0.1, add=T)
Я никогда не тратил так много времени, пытаясь эффективно векторизовать один единственный цикл в моей жизни! Кроме того, я не против создания пунтинга и просто делаю это с C++ и Rcpp, но я решил сначала спросить здесь, есть ли способ в R сделать его масштабируемым и быстрее. Заранее спасибо!
2 ответа
Вот решение, которое дает почти 100-кратное увеличение скорости. Я не до конца понимаю, почему улучшение настолько велико, но, возможно, один из реальных экспертов по таблицам данных может прокомментировать это.
library(RANN)
library(data.table)
set.seed(1L) # for reproducible data
DATA <- data.table(runif(10000, 0,1),
runif(10000, 0,1),
runif(10000, 0,1),
runif(10000, 10,30))
colnames(DATA)<-c("x","y","z","P")
nn.idx <- nn2(DATA[,1:3], DATA[,1:3], k=50,
treetype = "kd", searchtype = "radius",
radius = 0.075)$nn.idx
# (1)
# Timing for original loop.
system.time(for(index in 1:nrow(DATA)) {
DATA$mean.P[index] <- mean(DATA[nn.idx[index,], P])
})
# user system elapsed
# 7.830 0.850 8.684
# (2)
# Use `set()` instead of `$<-` and `[<-`.
system.time({for(index in 1:nrow(DATA)) {
set(DATA, i=index, j="mean_P_2", value=mean(DATA[nn.idx[index, ], P]))
}})
# user system elapsed
# 3.405 0.008 3.417
Как вы можете видеть, улучшение происходит в 2 раза, если подставить data.table-specific set() функция в оригинальном цикле.
Затем я попытался поместить всю функциональность в функции, специфичные для data.table (в основном внутри синтаксиса data.table []). Я также положил P значения в вектор, потому что доступ к значениям в векторах обычно намного быстрее, чем аналогичные операции над data.frames или data.tables.
# (3)
# Add row index.
DATA[, row_idx:=seq(nrow(DATA))]
# Isolate P values in a vector, because vector access is cheaper
# than data.table or data.frame access.
P_vec = DATA$P
system.time({
# Create a list column where each element is a vector of 50 integer indexes.
DATA[, nn_idx:=lapply(row_idx, function(i) nn.idx[i, ])]
# Use `:=` and `by=` to internalize the loop within `[.data.table`.
DATA[, mean_P_3:=mean(P_vec[nn_idx[[1]]]), by=row_idx]
})
# user system elapsed
# 0.092 0.002 0.095
# All results are identical.
all.equal(DATA$mean.P, DATA$mean_P_2)
# [1] TRUE
all.equal(DATA$mean.P, DATA$mean_P_3)
# [1] TRUE
Это дает почти 100-кратное увеличение скорости по сравнению с исходным циклом.
Похоже, масштабируется до 1 миллиона точек данных:
# Try with 1 million data points.
set.seed(1L) # for reproducible data
DATA2 <- data.table(runif(1e6, 0,1),
runif(1e6, 0,1),
runif(1e6, 0,1),
runif(1e6, 10,30))
colnames(DATA2) <- c("x","y","z","P")
system.time({
nn.idx2 <- nn2(DATA2[,1:3], DATA2[,1:3], k=50,
treetype = "kd", searchtype = "radius",
radius = 0.075)$nn.idx
})
# user system elapsed
# 346.603 1.883 349.708
DATA2[, row_idx:=seq(nrow(DATA2))]
P_vec = DATA2$P
system.time({
DATA2[, nn_idx:=lapply(row_idx, function(i) nn.idx2[i, ])]
DATA2[, mean_P:=mean(P_vec[nn_idx[[1]]]), by=row_idx]
})
# user system elapsed
# 15.685 0.587 16.297
Время было сделано на одном ядре MacBook Pro 2011 года (Sandy Bridge 2,2 ГГц). Использование оперативной памяти осталось ниже 1,5 ГБ.
Вот еще одно решение с использованием melt() чтобы изменить матрицу индекса в длинном формате, объединяя и агрегируя:
long <- melt(as.data.table(nn.idx)[, pt := .I], measure.vars = patterns("V"))
tmp <- long[DATA[, pt := .I], on = .(value = pt)][, mean(P), by = .(pt)][order(pt), V1]
DATA[, mean.P := tmp][, pt := NULL][]
объяснение
Индексная матрица nn.idx преобразуется в таблицу data.table и получает столбец pt который является идентификатором строки точек. Затем матрица преобразуется из широкого в длинный формат.
tmp вектор средних значений соседних точек. Они найдены путем правильного присоединения DATA с long сопоставить индексы ближайших соседних точек (в столбце value) с индексом точки, добавленным к DATA заранее.
Последний шаг - добавить результат в виде нового столбца в DATA,
Вариант 2
В качестве альтернативы, промежуточный результат может быть добавлен с помощью второго соединения:
long <- melt(as.data.table(nn.idx)[, pt := .I], measure.vars = patterns("V"))
long[DATA[, pt := .I], on = .(value = pt)][, mean(P), by = .(pt)][DATA, on = "pt"]