Как минимизировать неприемлемо длительное время выполнения созданного R-кода

ProgramBox

Как минимизировать неприемлемо длительное время выполнения созданного R-кода

Post author:admin
Запись опубликована:29 июня, 2022
Post category:Вопросы по программированию

#r #for-loop #runtime #missing-data #large-data

Вопрос:

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

Приведенная ниже структура на 100% соответствует тому, чего я пытаюсь достичь, когда я тестирую с очень небольшим количеством точек данных. Но поскольку количество столбцов и строк становится 2781 и 280 соответственно, я понимаю, что это занимает вечность, хотя я на 100% уверен, что это работает правильно, даже когда я вижу обновленное окно среды моей R-Studio каждый раз, когда я его обновляю.

В моих данных также много пропущенных значений, вероятно, 40% или что-то в этом роде. Я думаю, что это также увеличивает время вычислений.

Размер данных равен 315 * 2781. Тем не менее, я пытаюсь получить результат в форме матрицы 280 * 2781.

Могу ли я, пожалуйста, получить помощь в минимизации времени выполнения этого следующего кода? Было бы очень признательно, если я смогу!

 options(java.parameters = "- Xmx8000m")
memory.limit(size=8e 6)

data=read.table("C:/Data/input.txt",T,sep="t");
data=data.frame(data)[,-1]



corr<-NULL
corr2<-NULL
corr3<-NULL


for(i in 1:280) 
{
  corr2<-NULL
  for(j in 1:2781)
  {
    data2<-data[,-j]
    corr<-NULL
    for(k in 1:2780)
    {
      ifelse((is.error(grangertest(data[i:(i 35),j] ~ data2[i:(i 35),k], order = 1, na.action = na.omit)$P[2])==TRUE) || (grangertest(data[i:(i 35),j] ~ data2[i:(i 35),k], order = 1, na.action = na.omit)$P[2])>0.05|| (is.na(grangertest(data[i:(i 35),j] ~ data2[i:(i 35),k], order = 1, na.action = na.omit)$P[2])==TRUE),corr<-cbind(corr,0),corr<-cbind(corr,1))
    }
    corr2<-rbind(corr2,corr)
  }
  corr3<-rbind(corr3,rowSums(corr2))
}

Фрагмент моих данных выглядит следующим образом:

 > dput(data[1:30, 1:10])
structure(c(0.567388170165941, 0.193093325709924, 0.965938209090382, 
0.348295788047835, 0.496113050729036, 0.0645384560339153, 0.946750836912543, 
0.642093246569857, 0.565092500532046, 0.0952424583956599, 0.444063827162609, 
0.709971546428278, 0.756330407923087, 0.601746253203601, 0.341865634545684, 
0.953319212188944, 0.0788604547269642, 0.990508111426607, 0.35519331949763, 
0.697004508692771, 0.285368352662772, 0.274287624517456, 0.575733694015071, 
0.12937490013428, 0.00476219342090189, 0.684308280004188, 0.189448777819052, 
0.615732178557664, 0.404873769031838, 0.357331350911409, 0.565436001634225, 
0.380773033713922, 0.348490287549794, 0.0473814208526164, 0.389312234241515, 
0.562123290728778, 0.30642102798447, 0.911173274740577, 0.566258994862437, 
0.837928073247895, 0.107747194357216, 0.253737836843356, 0.651503744535148, 
0.187739939894527, 0.951192815322429, 0.740037888288498, 0.0817571650259197, 
0.740519099170342, 0.601534485351294, 0.120900869136676, 0.415282893227413, 
0.591146623482928, 0.698511375114322, 0.08557975362055, 0.139396222075447, 
0.303953414550051, 0.0743798329494894, 0.0293272000271827, 0.335832208395004, 
0.665010208031163, 0.0319741254206747, 0.678886031731963, 0.154593498911709, 
0.275712370406836, 0.828485634410754, 0.921500099124387, 0.651940459152684, 
0.00574865937232971, 0.82236105017364, 0.55089360428974, 0.209424041677266, 
0.861786168068647, 0.672873278381303, 0.301034058211371, 0.180336013436317, 
0.481560358777642, 0.901354183442891, 0.986482679378241, 0.90117057505995, 
0.476308439625427, 0.638073122361675, 0.27481731469743, 0.689271076582372, 
0.324349449947476, 0.56620552809909, 0.867861548438668, 0.78374840435572, 
0.0668482843320817, 0.276675389613956, 0.990600393852219, 0.990227151894942, 
0.417612489778548, 0.391012848122045, 0.348758921027184, 0.0799746725242585, 
0.88941288786009, 0.511429069796577, 0.0338982092216611, 0.240115304477513, 
0.0268365524243563, 0.67206134647131, 0.816803207853809, 0.344421110814437, 
0.864659120794386, 0.84128700569272, 0.116056860191748, 0.303730394458398, 
0.48192183743231, 0.341675494797528, 0.0622653553728014, 0.823110743425786, 
0.483212807681412, 0.968748248415068, 0.953057422768325, 0.116025703493506, 
0.327919023809955, 0.590675016632304, 0.832283023977652, 0.342327545629814, 
0.576901035616174, 0.942689201096073, 0.59300709143281, 0.565881528891623, 
0.600007816683501, 0.133237989619374, 0.873827134957537, 0.744597729761153, 
0.755133397178724, 0.0245723063126206, 0.97799762734212, 0.636845340020955, 
0.73828601022251, 0.644093665992841, 0.57204390084371, 0.496023115236312, 
0.703613247489557, 0.149237307952717, 0.0871439634356648, 0.0632112647872418, 
0.83703236351721, 0.433215840253979, 0.430483993608505, 0.924051651498303, 
0.913056606892496, 0.914889572421089, 0.215407102368772, 0.76880722376518, 
0.269207723205909, 0.865548757137731, 0.28798541566357, 0.391722843516618, 
0.649806497385725, 0.459413924254477, 0.907465039752424, 0.48731207777746, 
0.554472463205457, 0.779784266138449, 0.566323830280453, 0.208658932242543, 
0.958056638715789, 0.61858483706601, 0.838681482244283, 0.286310768220574, 
0.895410191034898, 0.448722236789763, 0.297688684659079, 0.33291415637359, 
0.0115265529602766, 0.850776052568108, 0.764857453294098, 0.469730701530352, 
0.222089925780892, 0.0496484278701246, 0.32886885642074, 0.356443469878286, 
0.612877089297399, 0.727906176587567, 0.0292073413729668, 0.429160050582141, 
0.232313714455813, 0.678631312213838, 0.642334033036605, 0.99107678886503, 
0.542449960019439, 0.835914565017447, 0.52798323193565, 0.303808332188055, 
0.919654499506578, 0.944237019168213, 0.52141259261407, 0.794379767496139, 
0.72268659202382, 0.114752230467275, 0.175116094760597, 0.437696389388293, 
0.852590200025588, 0.511136321350932, 0.30879021063447, 0.174206420546398, 
0.14262041519396, 0.375411552377045, 0.0204910831525922, 0.852757754037157, 
0.631567053496838, 0.475924106314778, 0.508682047016919, 0.307679089019075, 
0.70284536993131, 0.851252349093556, 0.0868967010173947, 0.586291917832568, 
0.0529140203725547, 0.440692059928551, 0.207642213441432, 0.777513341512531, 
0.141496006632224, 0.548626560717821, 0.419565241318196, 0.0702310993801802, 
0.499403427587822, 0.189343606121838, 0.370725362794474, 0.888076487928629, 
0.83070912421681, 0.466137421084568, 0.177098380634561, 0.91202046489343, 
0.142300580162555, 0.823691181838512, 0.41561916610226, 0.939948018174618, 
0.806491429451853, 0.795849160756916, 0.566376683535054, 0.36814984655939, 
0.307756055146456, 0.602875682059675, 0.506007500691339, 0.538658684119582, 
0.420845189364627, 0.663071365095675, 0.958144341595471, 0.793743418296799, 
0.983086514985189, 0.266262857476249, 0.817585011478513, 0.122843299992383, 
0.989197303075343, 0.71584410732612, 0.500571243464947, 0.397394519997761, 
0.659465527161956, 0.459530522814021, 0.602246116613969, 0.250076721422374, 
0.17533828667365, 0.6599256307818, 0.184704560553655, 0.15679649473168, 
0.513444944983348, 0.205572377191857, 0.430164282443002, 0.131548407254741, 
0.914019819349051, 0.935795902274549, 0.857401241315529, 0.977940042736009, 
0.41389597626403, 0.179183913161978, 0.431347143370658, 0.477178965462372, 
0.121315707685426, 0.107695729471743, 0.634954946814105, 0.859707030234858, 
0.855825762730092, 0.708672808250412, 0.674073817208409, 0.672288877889514, 
0.622144045541063, 0.433355041313916, 0.952878215815872, 0.229569894727319, 
0.289388840552419, 0.937473804224283, 0.116283216979355, 0.659604362910613, 
0.240837284363806, 0.726138337515295, 0.68390148691833, 0.381577257299796, 
0.899390475358814, 0.26472729514353, 0.0383855854161084, 0.855232689995319, 
0.655799814499915, 0.335587574867532, 0.163842789363116, 0.0353666560258716, 
0.048316186061129), .Dim = c(30L, 10L))

1. Было бы здорово прокомментировать код по пути, чтобы убедиться, что мы понимаем, что он должен делать. В качестве общего предложения в R: — c , rbind , cbind переписывать векторы очень дорого с точки зрения вычислений, лучше всего было бы создать выходной вектор с ожидаемым количеством элементов, а затем перезаписать его в каждой отдельной позиции с помощью выходных данных цикла for . — поскольку циклы в R очень медленные, я предлагаю вам использовать apply sapply lapply синтаксис // — удаление одной строки за раз ( data[,-j] ) также может быть довольно тяжелым — удаление одного столбца за раз также

2. Если на вашей машине несколько ядер, вы можете распараллелить свой код, @Eric. См.: Распараллеленные циклы с R

3. Не data[i:(i 35),j] будет ли ошибка после i получения 246 ?

4. Можете ли вы опубликовать образцы данных? Пожалуйста, отредактируйте вопрос с выводом dput(data[1:30, 1:10]) .

5. @RuiBarradas: Пожалуйста, проверьте мой dput выше.

Ответ №1:

Я преобразовал только внутренний цикл в mapply и провел быстрый тест скорости:

 library(lmtest)

data <- matrix(runif(315*2781), nrow = 315)

get01 <- function(x, y) {
  try(gt <- grangertest(x ~ y, order = 1, na.action = na.omit)$P[2])
  if (exists("gt")) {
    if (gt > 0.05 || is.na(gt)) {
      return(0)
    } else {
      return(1)
    }
  } else {
    return(0)
  }
}

i <- 1; j <- 1
system.time(corr <- mapply(function(k) {get01(data[i:(i 35),j], data[i:(i 35),k])}, (1:2781)[-j]))
#>    user  system elapsed 
#>  21.505   0.014  21.520

Для этого потребуется выполнить это mapply 778680 раз, так что это займет около 200 дней. Вам понадобится либо другой подход с тестом Грейнджера, либо несколько ядер. Вот команда для замены полного цикла:

 corr3 <- t(mapply(function(i) colSums(mapply(function(j) mapply(function(k) {get01(data[i:(i 35),j], data[i:(i 35),k])}, (1:2781)[-j]), 1:2781)), 1:280))

Замените это сначала mapply на simplify2array(parLapply для распараллеливания:

 library(parallel)

cl <- makeCluster(detectCores())
clusterExport(cl, list("data", "get01"))
parLapply(cl, cl, function(x) require(lmtest))
corr3 <- t(simplify2array(parLapply(cl, 1:280, function(i) colSums(mapply(function(j) mapply(function(k) {get01(data[i:(i 35),j], data[i:(i 35),k])}, (1:2781)[-j]), 1:2781)))))
stopCluster(cl)

1. Я перейду на компьютер с Windows и разработаю parLappply. Однако, чтобы добраться до него, может потребоваться некоторое время. Спасибо за разъяснение размера строки.

2. parLapply не так сложно использовать. Я запустил свой компьютер с Windows и тестирую код.

3. Просто любопытно — сколько у вас ядер (из detectCores() )?

4. Я обновил этот пост с parLapply помощью . Если Rui granger_test работает для вас, использование этого с тем, что у меня есть выше, вероятно, будет близко к лучшему решению без полного пересмотра lmtest::grangertest , и lmtest::waldtest это позволит вам векторизовать вашу матрицу.

5. Просто обновил его, должно было быть parLapply(cl, cl, function(x) require(lmtest))

Ответ №2:

Вот версия, не распараллеленная, которая ускоряет код в вопросе в разы больше, чем в 4.

Некоторые узкие места в коде вопроса легко обнаружить:

Матрицы corr? расширяются внутри циклов. Решение состоит в том, чтобы заранее зарезервировать память;
Тест grangertest вызывается 3 раза за внутреннюю итерацию, когда требуется только один;
Для cbind с 0 или 1 фактически создает вектор, а не матрицу.

Вот сравнительный тест между кодом вопроса и приведенной ниже функцией.

 library(lmtest)

# avoids loading an extra package
is.error <- function(x){
  inherits(x, c("error", "try-error"))
}

Lag <- 5L
nr <- nrow(data)
nc <- ncol(data)

t0 <- system.time({
  corr<-NULL
  corr2<-NULL
  corr3<-NULL
  for(i in 1:(nr - Lag))
  {
    corr2<-NULL
    data3 <- data[i:(i   Lag), ]
    for(j in 1:nc)
    {
      data2<-data[,-j]
      corr<-NULL
      for(k in 1:(nc - 1L))
      {
        ifelse((is.error(grangertest(data[i:(i Lag),j] ~ data2[i:(i Lag),k], order = 1, na.action = na.omit)$P[2])==TRUE) ||
                 (grangertest(data[i:(i Lag),j] ~ data2[i:(i Lag),k], order = 1, na.action = na.omit)$P[2])>0.05 ||
                 (is.na(grangertest(data[i:(i Lag),j] ~ data2[i:(i Lag),k], order = 1, na.action = na.omit)$P[2])==TRUE),
               corr<-cbind(corr,0),
               corr<-cbind(corr,1)
        )
      }
      corr2 <- rbind(corr2, corr)
    }
    corr3<-rbind(corr3, rowSums(corr2))
  }
  corr3
})

Я буду использовать упрощенную версию lmtest::grangertest .

 granger_test <- function (x, y, order = 1, na.action = na.omit, ...) {
  xnam <- deparse(substitute(x))
  ynam <- deparse(substitute(y))
  n <- length(x)
  all <- cbind(x = x[-1], y = y[-1], x_1 = x[-n], y_1 = y[-n])
  y <- as.vector(all[, 2])
  lagX <- as.matrix(all[, (1:order   2)])
  lagY <- as.matrix(all[, (1:order   2   order)])
  fm <- lm(y ~ lagY   lagX)
  rval <- lmtest::waldtest(fm, 2, ...)
  attr(rval, "heading") <- c("Granger causality testn", paste("Model 1: ", 
                                                               ynam, " ~ ", "Lags(", ynam, ", 1:", order, ")   Lags(", 
                                                               xnam, ", 1:", order, ")nModel 2: ", ynam, " ~ ", "Lags(", 
                                                               ynam, ", 1:", order, ")", sep = ""))
  rval
}

А теперь функция для запуска тестов.

 f_Rui <- function(data, Lag){
  nr <- nrow(data)
  nc <- ncol(data)
  corr3 <- matrix(0, nrow = nr - Lag, ncol = nc)
  data3 <- matrix(0, nrow = Lag   1L, ncol = nc)
  data2 <- matrix(0, nrow = Lag   1L, ncol = nc - 1L)
  for(i in 1:(nr - Lag)) {
    corr2 <- matrix(0, nrow = nc, ncol = nc - 1L)
    data3[] <- data[i:(i   Lag), ]
    for(j in 1:nc) {
      corr <- integer(nc - 1L)
      data2[] <- data3[, -j]
      for(k in 1:(nc - 1L)){
        res <- tryCatch(
          grangertest(x = data2[, k], y = data3[, j], order = 1, na.action = na.omit),
          error = function(e) e
        )
        if(!inherits(res, "error") amp;amp; !is.na(res[['Pr(>F)']][2]) amp;amp; res[['Pr(>F)']][2] <= 0.05) {
          corr[k] <- 1L
        }
      }
      corr2[j, ] <- corr
    }
    corr3[i, ] <- rowSums(corr2)
  }
  corr3
}

Результаты identical и сроки намного лучше.

 t1 <- system.time({
  res <- f_Rui(data, 5L)
})

identical(corr3, res)
#[1] TRUE

times <- rbind(t0, t1)
t(t(times)/t1)
#   user.self sys.self  elapsed user.child sys.child
#t0  4.682908 1.736111 4.707783        NaN       NaN
#t1  1.000000 1.000000 1.000000        NaN       NaN

1. Баррадс: Спасибо, но я думаю, что «f_Rui (data, 5L)» на последнем шаге, к сожалению, не выполняется для моего случая. В нем говорится: «Ошибка в f_Rui (данные, 5L): неиспользуемый аргумент (5)» Кроме того, не могли бы вы объяснить, что такое 5L и почему он задан? Большое спасибо.

2. Извините, я вижу, что задержка — это то, что я настроил как скользящее окно, как мне понравилось 35.

3. Спасибо! Это была хорошая возможность учиться!

Метки: Как свести к минимуму недопустимо длительное время выполнения созданного кода R