Перемещение за пределы функции R optim

Я пытаюсь использовать R для оценки многочленной логитной модели с ручной спецификацией. Я нашел несколько пакетов, которые позволяют вам оценивать модели MNL здесь или здесь.

Я нашел несколько других статей о "прокатке" вашей собственной функции MLE здесь. Однако, из-за моего копания - все эти функции и пакеты полагаются на внутреннюю функцию optim.

В моих тестовых тестах optim является узким местом. Используя смоделированный набор данных с ~ 16000 наблюдениями и 7 параметрами, R занимает около 90 секунд на моей машине. Эквивалентная модель в Biogeme занимает ~ 10 секунд. Коллега, который пишет свой собственный код в Ox, сообщает около 4 секунд для этой же модели.

Есть ли у кого-нибудь опыт написания собственной функции MLE или можно указать мне в сторону чего-то, что оптимизировано за пределами функции optim по умолчанию (каламбур не предназначен)?

Если кто-то хочет, чтобы R-код воссоздал модель, дайте мне знать - я получу ее. Я не предоставил его, поскольку он не имеет прямого отношения к проблеме оптимизации функции optim и для сохранения пространства...

EDIT: Спасибо всем за ваши мысли. Основываясь на множестве комментариев ниже, мы смогли получить R в том же самом поле, что и Biogeme, для более сложных моделей, и R действительно быстрее для нескольких более мелких/более простых моделей, которые мы запускали. Я думаю, что долгосрочное решение этой проблемы будет связано с написанием отдельной функции максимизации, которая опирается на fortran или библиотеку C, но, безусловно, открыта для других подходов.

Ответ 1

Пробовал с функцией nlm() уже? Не знаю, если это намного быстрее, но оно улучшает скорость. Также проверьте параметры. Оптимист использует медленный алгоритм по умолчанию. Вы можете получить > 5-кратное ускорение, используя алгоритм Квази-Ньютона (метод = "BFGS" ) вместо значения по умолчанию. Если вы не слишком беспокоитесь о последних цифрах, вы также можете установить уровни допуска выше nlm(), чтобы получить дополнительную скорость.

f <- function(x) sum((x-1:length(x))^2)

a <- 1:5

system.time(replicate(500,
     optim(a,f)
))
   user  system elapsed 
   0.78    0.00    0.79 

system.time(replicate(500,
     optim(a,f,method="BFGS")
))
   user  system elapsed 
   0.11    0.00    0.11 

system.time(replicate(500,
     nlm(f,a)
))
   user  system elapsed 
   0.10    0.00    0.09 

system.time(replicate(500,
      nlm(f,a,steptol=1e-4,gradtol=1e-4)
))
   user  system elapsed 
   0.03    0.00    0.03

Ответ 2

Вы рассмотрели материал в Просмотр задачи CRAN для оптимизации?

Ответ 3

FWIW, я сделал это в C-ish, используя OPTIF9. Тебе было бы тяжело идти быстрее. Существует множество способов замедления работы, например, путем запуска интерпретатора, такого как R.

Добавлено: Из комментариев видно, что OPTIF9 используется в качестве механизма оптимизации. Это означает, что, скорее всего, основная часть времени тратится на оценку целевой функции в R. Хотя возможно, что функции C используются под некоторыми операциями, все же накладные расходы переводчика. Существует быстрый способ определить, какие строки вызовов кода и функций в R ответственны в течение большей части времени, а именно, чтобы приостановить его с помощью клавиши Escape и изучить стек. Если оператор стоит X% времени, он находится в стеке X% времени. Вы можете обнаружить, что есть операции, которые не идут на C и должны быть. Любой коэффициент ускорения, который вы получите таким образом, будет сохранен, когда вы найдете способ распараллеливать выполнение R.

Ответ 4

Я являюсь автором пакета R optimParallel, который может быть полезен в вашем случае. Пакет предоставляет параллельные версии методов оптимизации на основе градиента optim(). Основной функцией пакета является optimParallel(), который имеет одинаковое использование и вывод, как optim(). Использование optimParallel() может значительно сократить время оптимизации, как показано на следующем рисунке (p - количество параметров). Для получения дополнительной информации см. Https://cran.r-project.org/package=optimParallel и http://arxiv.org/abs/1804.11058.