Bigglm в R: ограничения и улучшения исходного кода (например, вызов Fortran)

Недавно я использовал R для запуска обобщенной линейной модели (GLM) в файле csv размером 100 МБ (9 миллионов строк по 5 столбцов). Содержимое этого файла включает 5 столбцов, называемых depvar, var1, var2, var3, var4, и все они распределены случайным образом, так что столбцы содержат числа 0,1 или 2. В основном я использовал пакет biglm для запуска GLM на этот файл данных, и R обработал его примерно за 2 минуты. Это было на Linux-машине с R версии 2.10 (сейчас я обновляюсь до 2.14), 4 ядрами и 8 ГБ ОЗУ. Обычно я хочу запускать код быстрее, примерно за 30-60 секунд. Одним из решений является добавление большего количества ядер и оперативной памяти, но это будет только временное решение, поскольку я понимаю, что наборы данных будут только увеличиваться. В идеале я хочу найти способ ускорить код для bigglm. Я запустил некоторый код профилирования R. Добавляем следующий код (перед кодом, который я хочу запустить, чтобы проверить его скорость):

Rprof('method1.out')

Затем, набрав эту команду, я пишу свой код bigglm, который выглядит примерно так:

x<-read.csv('file location of 100mb file')
form<-depvar~var1+var2+var3+var4
a<-bigglm(form, data=x, chunksize=800000, sandwich=FALSE, family=binomial())
summary(a)
AIC(a)
deviance(a)

После запуска этих кодов, которые занимают от 2 до 3 минут, я набираю следующее, чтобы увидеть свой код профилирования:

Rprofsummary('method1.out')

Затем я получаю разбивку процесса bigglm и то, какие отдельные строки занимают очень много времени. После просмотра я был удивлен, увидев, что был вызов кода fortran, который занимал очень много времени (около 20 секунд). Этот код можно найти в базовом файле Bigglm по адресу:

http://cran.r-project.org/web/packages/biglm/index.html

В файле bigglm 0.8.tar.gz

В основном я спрашиваю сообщество, можно ли сделать этот код быстрее? Например, изменив код, чтобы вызвать код Fortran для выполнения QR-разложения. Кроме того, были другие функции, такие как as.character и model.matrix, которые также занимали много времени. Я не прилагал сюда файл профилирования, так как считаю, что его можно легко воспроизвести, учитывая предоставленную мной информацию, но в основном я намекаю на большую проблему больших данных и обработки GLM на этих больших данных. Это проблема, которую разделяет сообщество R, и я думаю, что любые отзывы или помощь будут благодарны по этому вопросу. Вероятно, вы можете легко воспроизвести этот пример, используя другой набор данных, и посмотреть, что занимает так много времени в коде bigglm, и увидеть, совпадают ли они с тем, что я нашел. Если да, может кто-нибудь, пожалуйста, помогите мне выяснить, как заставить bigglm работать быстрее. После того, как Бен запросил это, я загрузил фрагмент кода профилирования, который у меня был, а также первые 10 строк моего файла csv:

 CSV File:

 var1,var2,var3,var4,depvar
   1,0,2,2,1
   0,0,1,2,0
   1,1,0,0,1
   0,1,1,2,0
   1,0,0,3,0
   0,0,2,2,0
   1,1,0,0,1
   0,1,2,2,0
   0,0,2,2,1

Этот вывод CSV был скопирован из моего текстового редактора UltraEdit, и видно, что var1 принимает значения 0 или 1, var2 принимает значения 0 и 1, var3 принимает значения 0,1,2, var4 принимает значения 0,1, 2,3 и depvar принимает значения 1 или 0. Этот CSV может быть воспроизведен в Excel с помощью функции RAND примерно до 1 миллиона строк, затем его можно скопировать и вставить несколько раз, чтобы получить большое количество строк в текстовом редакторе. как ultraedit. Обычно введите RAND () в один столбец для 1 миллиона столбцов, затем выполните округление (столбец) в столбце рядом со столбцом RAND (), чтобы получить единицы и нули. То же самое относится и к 0,1,2,3.

Файл профилирования длинный, поэтому я приложил строки, которые занимали больше всего времени:

 summaryRprof('method1.out')



$by.self
                        self.time self.pct total.time total.pct
"model.matrix.default"      25.40     20.5      26.76      21.6
".Call"                     20.24     16.3      20.24      16.3
"as.character"              17.22     13.9      17.22      13.9
"[.data.frame"              14.80     11.9      22.12      17.8
"update.bigqr"               5.72      4.6      14.38      11.6
"-"                          4.36      3.5       4.36       3.5
"anyDuplicated.default"      4.18      3.4       4.18       3.4
"|"                          3.98      3.2       3.98       3.2
"*"                          3.44      2.8       3.44       2.8
"/"                          3.18      2.6       3.18       2.6
"unclass"                    2.28      1.8       2.28       1.8
"sum"                        2.26      1.8       2.26       1.8
"attr"                       2.12      1.7       2.12       1.7
"na.omit"                    2.02      1.6      20.00      16.1
"%*%"                        1.74      1.4       1.74       1.4
"^"                          1.56      1.3       1.56       1.3
"bigglm.function"            1.42      1.1     122.52      98.8
"+"                          1.30      1.0       1.30       1.0
"is.na"                      1.28      1.0       1.28       1.0
"model.frame.default"        1.20      1.0      22.68      18.3
">"                          0.84      0.7       0.84       0.7
"strsplit"                   0.62      0.5       0.62       0.5

$by.total
                        total.time total.pct self.time self.pct
"standardGeneric"           122.54      98.8      0.00      0.0
"bigglm.function"           122.52      98.8      1.42      1.1
"bigglm"                    122.52      98.8      0.00      0.0
"bigglm.data.frame"         122.52      98.8      0.00      0.0
"model.matrix.default"       26.76      21.6     25.40     20.5
"model.matrix"               26.76      21.6      0.00      0.0
"model.frame.default"        22.68      18.3      1.20      1.0
"model.frame"                22.68      18.3      0.00      0.0
"["                          22.14      17.9      0.02      0.0
"[.data.frame"               22.12      17.8     14.80     11.9
".Call"                      20.24      16.3     20.24     16.3
"na.omit"                    20.00      16.1      2.02      1.6
"na.omit.data.frame"         17.98      14.5      0.02      0.0
"model.response"             17.44      14.1      0.10      0.1
"as.character"               17.22      13.9     17.22     13.9
"names<-"                    17.22      13.9      0.00      0.0
"<Anonymous>"                15.10      12.2      0.00      0.0
"update.bigqr"               14.38      11.6      5.72      4.6
"update"                     14.38      11.6      0.00      0.0
"data"                       10.26       8.3      0.00      0.0
"-"                           4.36       3.5      4.36      3.5
"anyDuplicated.default"       4.18       3.4      4.18      3.4
"anyDuplicated"               4.18       3.4      0.00      0.0
"|"                           3.98       3.2      3.98      3.2
"*"                           3.44       2.8      3.44      2.8
"/"                           3.18       2.6      3.18      2.6
"lapply"                      3.04       2.5      0.04      0.0
"sapply"                      2.44       2.0      0.00      0.0
"as.list.data.frame"          2.30       1.9      0.02      0.0
"as.list"                     2.30       1.9      0.00      0.0
"unclass"                     2.28       1.8      2.28      1.8
"sum"                         2.26       1.8      2.26      1.8
"attr"                        2.12       1.7      2.12      1.7
"etafun"                      1.88       1.5      0.14      0.1
"%*%"                         1.74       1.4      1.74      1.4
"^"                           1.56       1.3      1.56      1.3
"summaryRprof"                1.48       1.2      0.02      0.0
"+"                           1.30       1.0      1.30      1.0
"is.na"                       1.28       1.0      1.28      1.0
">"                           0.84       0.7      0.84      0.7
"FUN"                         0.70       0.6      0.36      0.3
"strsplit"                    0.62       0.5      0.62      0.5

Больше всего меня удивила функция .Call, которая вызывает Фортран. Может, я этого не понял. Похоже, что после использования этой функции все вычисления производятся. Я думал, что это похоже на функцию связывания, которая использовалась для извлечения кода Fortran. Кроме того, если Фортран выполняет всю работу и все итеративно взвешенные наименьшие квадраты / QR, почему тогда остальная часть кода занимает так много времени.