Ускорение интерполяции
Я провожу около 45 000 локальных линейных регрессий (по сути) около 1,2 миллиона наблюдений, поэтому я был бы признателен за помощь в попытке ускорить процесс, потому что я нетерпелив.
Я в основном пытаюсь составить контракты по заработной плате за год - функция заработной платы (опыт с учетом фирмы, год, должность)- для нескольких фирм.
Вот набор данных (базовая структура), с которым я работаю:
> wages
firm year position exp salary
1: 0007 1996 4 1 20029
2: 0007 1996 4 1 23502
3: 0007 1996 4 1 22105
4: 0007 1996 4 2 23124
5: 0007 1996 4 2 22700
---
1175141: 994 2012 5 2 47098
1175142: 994 2012 5 2 45488
1175143: 994 2012 5 2 47098
1175144: 994 2012 5 3 45488
1175145: 994 2012 5 3 47098
Я хочу построить функцию заработной платы для уровней опыта от 0 до 40 для всех фирм, а-ля:
> salary_scales
firm year position exp salary
1: 0007 1996 4 0 NA
2: 0007 1996 4 1 21878.67
3: 0007 1996 4 2 23401.33
4: 0007 1996 4 3 23705.00
5: 0007 1996 4 4 24260.00
---
611019: 9911 2015 4 36 NA
611020: 9911 2015 4 37 NA
611021: 9911 2015 4 38 NA
611022: 9911 2015 4 39 NA
611023: 9911 2015 4 40 NA
С этой целью я работал (по предложению @BondedDust здесь) с пакетом COBS (COnstrained B-Spline), который позволяет мне встроить монотонность контракта о заработной плате.
Некоторые проблемы остаются; в частности, когда мне нужно экстраполировать (когда в данной фирме нет очень молодых или очень старых сотрудников), наблюдается тенденция к потере монотонности или падению ниже 0.
Чтобы обойти это, я использовал простую линейную экстраполяцию вне границ данных - расширить кривую подгонки за пределы min_exp а также max_exp так что он проходит через две самые низкие (или самые высокие) точки подгонки - не идеально, но, похоже, все идет хорошо.
Имея это в виду, вот как я делаю это до сих пор (имейте в виду, что я data.table фанатик):
#get the range of experience for each firm
wages[,min_exp:=min(exp),by=.(year,firm,position)]
wages[,max_exp:=max(exp),by=.(year,firm,position)]
#Can't interpolate if there are only 2 or 3 unique experience cells represented
wages[,node_count:=length(unique(exp)),by=.(year,firm,position)]
#Nor if there are too few teachers
wages[,ind_count:=.N,by=.(year,firm,position)]
#Also troublesome when there is little variation in salaries like so:
wages[,sal_scale_flag:=mean(abs(salary-mean(salary)))<50,by=.(year,firm,position)]
wages[,sal_count_flag:=length(unique(salary))<5,by=.(year,firm,position)]
cobs_extrap<-function(exp,salary,min_exp,max_exp,
constraint="increase",print.mesg=F,nknots=8,
keep.data=F,maxiter=150){
#these are passed as vectors
min_exp<-min_exp[1]
max_exp<-min(max_exp[1],40)
#get in-sample fit
in_sample<-predict(cobs(x=exp,y=salary,
constraint=constraint,
print.mesg=print.mesg,nknots=nknots,
keep.data=keep.data,maxiter=maxiter),
z=min_exp:max_exp)[,"fit"]
#append by linear extension below min_exp
c(if (min_exp==1) NULL else in_sample[1]-
(min_exp:1)*(in_sample[2]-in_sample[1]),in_sample,
#append by linear extension above max_exp
if (max_exp==40) NULL else in_sample[length(in_sample)]+(1:(40-max_exp))*
(in_sample[length(in_sample)]-in_sample[length(in_sample)-1]))
}
salary_scales<-
wages[node_count>=7&ind_count>=10
&sal_scale_flag==0&sal_count_flag==0,
.(exp=0:40,
salary=cobs_extrap(exp,salary,min_exp,max_exp)),
by=.(year,firm,position)]
Заметили что-нибудь конкретное, что могло бы замедлить мой код? Или я вынужден быть терпеливым?
Чтобы поиграть, вот некоторые из небольших комбо-позиций:
firm year position exp salary count
1: 0063 2010 5 2 37433 10
2: 0063 2010 5 2 38749 10
3: 0063 2010 5 4 38749 10
4: 0063 2010 5 8 42700 10
5: 0063 2010 5 11 47967 10
6: 0063 2010 5 15 50637 10
7: 0063 2010 5 19 51529 10
8: 0063 2010 5 23 50637 10
9: 0063 2010 5 33 52426 10
10: 0063 2010 5 37 52426 10
11: 9908 2006 4 1 26750 10
12: 9908 2006 4 6 36043 10
13: 9908 2006 4 7 20513 10
14: 9908 2006 4 8 45023 10
15: 9908 2006 4 13 33588 10
16: 9908 2006 4 15 46011 10
17: 9908 2006 4 15 37179 10
18: 9908 2006 4 22 43704 10
19: 9908 2006 4 28 56078 10
20: 9908 2006 4 29 44866 10
1 ответ
В вашем коде есть много вещей, которые можно улучшить, но давайте сосредоточимся здесь на основном узком месте. Рассматриваемая проблема может рассматриваться как смущающая параллельная проблема. Это означает, что ваши данные могут быть разбиты на несколько более мелких частей, каждый из которых может быть вычислен в отдельных потоках по отдельности без каких-либо дополнительных затрат.
Чтобы увидеть возможности распараллеливания для текущей проблемы, вы должны сначала заметить, что вы выполняете точно такие же вычисления для каждой из отдельных фирм и / или лет в отдельности. Например, вы можете разделить вычисления на более мелкие подзадачи для каждого отдельного года, а затем распределить эти подзадачи по разным ядрам CPU/GPU. Таким образом, можно получить значительный выигрыш в производительности. Наконец, когда обработка подзадач завершена, единственное, что вам еще нужно сделать, это объединить результаты.
Однако R и все его внутренние библиотеки работают как единый поток. Вам придется явно разделить ваши данные и затем назначить подзадачи различным ядрам. Для достижения этого существует несколько пакетов R, которые поддерживают многопоточность. Мы будем использовать doparallel пакет в нашем примере здесь.
Вы не предоставили явный набор данных, достаточно большой для эффективного тестирования производительности, поэтому мы сначала создадим несколько случайных данных:
set.seed(42)
wages<-data.table(firm=substr(10001:10010,2,5)[sample(10,size=1e6,replace=T)],
year=round(unif(1e6,1996,2015)),
position=round(runif(1e6,4,5)),
exp=round(runif(1e6,1,40)),
salary=round(exp(rnorm(1e6,mean=10.682,sd=.286))))
> wages
firm year position exp salary
1: 0001 1996 4 14 66136
2: 0001 1996 4 3 42123
3: 0001 1996 4 9 46528
4: 0001 1996 4 11 35195
5: 0001 1996 4 2 43926
---
999996: 0010 2015 5 11 43140
999997: 0010 2015 5 23 64025
999998: 0010 2015 5 31 35266
999999: 0010 2015 5 11 36267
1000000: 0010 2015 5 7 44315
Теперь давайте запустим первую часть вашего кода:
#get the range of experience for each firm
wages[,min_exp:=min(exp),by=.(year,firm,position)]
wages[,max_exp:=max(exp),by=.(year,firm,position)]
#Can't interpolate if there are only 2 or 3 unique experience cells represented
wages[,node_count:=length(unique(exp)),by=.(year,firm,position)]
#Nor if there are too few teachers
wages[,ind_count:=.N,by=.(year,firm,position)]
#Also troublesome when there is little variation in salaries like so:
wages[,sal_scale_flag:=mean(abs(salary-mean(salary)))<50,by=.(year,firm,position)]
wages[,sal_count_flag:=length(unique(salary))<5,by=.(year,firm,position)]
> wages
firm year position exp salary min_exp max_exp node_count ind_count sal_scale_flag sal_count_flag
1: 0001 1996 4 14 66136 1 40 40 1373 FALSE FALSE
2: 0001 1996 4 3 42123 1 40 40 1373 FALSE FALSE
3: 0001 1996 4 9 46528 1 40 40 1373 FALSE FALSE
4: 0001 1996 4 11 35195 1 40 40 1373 FALSE FALSE
5: 0001 1996 4 2 43926 1 40 40 1373 FALSE FALSE
---
999996: 0010 2015 5 11 43140 1 40 40 1326 FALSE FALSE
999997: 0010 2015 5 23 64025 1 40 40 1326 FALSE FALSE
999998: 0010 2015 5 31 35266 1 40 40 1326 FALSE FALSE
999999: 0010 2015 5 11 36267 1 40 40 1326 FALSE FALSE
1000000: 0010 2015 5 7 44315 1 40 40 1326 FALSE FALSE
Теперь мы обработаем wages однопоточным способом, как вы делали раньше. Обратите внимание, что сначала мы сохраняем исходные данные, чтобы потом можно было выполнять многопоточные операции с ними и сравнивать результаты:
start <- Sys.time()
salary_scales_1 <-
wages[node_count>=7&ind_count>=10
&sal_scale_flag==0&sal_count_flag==0,
.(exp=0:40,salary=cobs_extrap(exp,salary,min_exp,max_exp)),
by=.(firm,year,position)]
print(paste("No Parallelisation time: ",Sys.time()-start))
> print(paste("No Parallelisation time: ",Sys.time()-start))
[1] "No Parallelisation time: 1.13971961339315"
> salary_scales_1
firm year position exp salary
1: 0001 1996 4 0 43670.14
2: 0001 1996 4 1 43674.00
3: 0001 1996 4 2 43677.76
4: 0001 1996 4 3 43681.43
5: 0001 1996 4 4 43684.99
---
16396: 0010 2015 5 36 44464.02
16397: 0010 2015 5 37 44468.60
16398: 0010 2015 5 38 44471.35
16399: 0010 2015 5 39 44472.27
16400: 0010 2015 5 40 43077.70
На обработку всего ушло около 1 минуты 8 секунд. Обратите внимание, что в нашем фиктивном примере у нас всего 10 разных фирм, поэтому время обработки не так существенно по сравнению с вашими локальными результатами.
Теперь давайте попробуем выполнить эту задачу параллельно. Как уже упоминалось, для нашего примера мы хотели бы разделить данные за год и назначить меньшие части отдельным ядрам. Мы будем использовать doParallel Пакет для этой цели:
Первое, что нам нужно сделать, это создать кластер с определенным количеством ядер. В нашем примере мы попытаемся использовать все доступные ядра. Далее нам нужно зарегистрировать кластер и экспортировать некоторые переменные в глобальные среды подузлов. В этом случае подузлам нужен только доступ к wages, Кроме того, некоторые из зависимых библиотек также должны быть оценены на узлах, чтобы заставить его работать. В этом случае узлам необходим доступ к обоим data.frame а также cobs библиотеки. Код выглядит так:
library(doParallel)
start <- Sys.time()
cl <- makeCluster(detectCores());
registerDoParallel(cl);
clusterExport(cl,c("wages"),envir=environment());
clusterEvalQ(cl,library("data.table"));
clusterEvalQ(cl,library("cobs"));
salary_scales_2 <- foreach(i = 1996:2015) %dopar%
{
subSet <- wages[.(i)] # binary subsetting
subSet[node_count>=7&ind_count>=10
&sal_scale_flag==0&sal_count_flag==0,
.(exp=0:40,
salary=cobs_extrap(exp,salary,min_exp,max_exp)),
by=.(firm,year,position)]
}
stopCluster(cl)
print(paste("With parallelisation time: ",Sys.time()-start))
> print(paste("With parallelisation time: ",Sys.time()-start))
[1] "With parallelisation time: 23.4177722930908"
Теперь у нас есть список таблиц данных salary_scales_2 который содержит подрезультаты для каждого отдельного года. Обратите внимание на ускорение обработки: на этот раз потребовалось всего 23 секунды вместо первоначальных 1,1 минуты (улучшение на 65%). Единственное, что нам еще нужно сделать сейчас - это объединить результаты. Мы можем использовать do.call("rbind", salary_scales_2) для того, чтобы объединить строки таблиц вместе (это почти не занимает времени - 002 секунды за один прогон). Наконец, мы также выполняем небольшую проверку, чтобы убедиться, что многопоточные результаты действительно идентичны результатам однопоточного прогона:
salary_scales_2<-do.call("rbind",salary_scales_2)
identical(salary_scales_1,salary_scales_2)
> identical(salary_scales_1,salary_scales_2)
[1] TRUE
ОТВЕТИТЬ НА КОММЕНТАРИЙ Это действительно очень интересный пример, но я думаю, что вы можете пропустить более важный вопрос здесь. data.table действительно выполняет оптимизацию, связанную с памятью и структурой, чтобы вы могли более эффективно запрашивать и получать доступ к вашим данным. Однако в этом примере отсутствует узкое место в основной памяти или поиске, особенно если вы сравниваете его с фактическим общим временем обработки данных в cobs функция. Например, строка, которую вы изменили subSet <- wages[year==uniqueYears[i],] занимает всего около 0,04 секунды на звонок, когда вы его время.
Если вы используете профилировщик на своих трассах, то вы заметите, что это не data.table или любой из его операций или группировок, которые просят о распараллеливании, это cobs функция, которая занимает почти все время обработки (и эта функция даже не использует data.table в качестве входных данных). В этом примере мы пытаемся переназначить нашу общую нагрузку на cobs функционировать на разных ядрах, чтобы добиться нашего ускорения. Наше намерение не состоит в том, чтобы разделить data.table операции, так как они совсем не дороги. Тем не менее, мы действительно должны разделить нашу data.table в результате того, что нам нужно разделить данные для отдельных cobs функция работает. На самом деле, мы даже воспользовались тем, что data.table эффективен во всех отношениях при разделении и объединении таблиц. Это не заняло дополнительного времени вообще.