Über die optim-Funktion von R hinausgehen

Question

Ich versuche, mit R ein multinomiales Logit-Modell mit einer manuellen Spezifikation zu schätzen. Ich habe ein paar Pakete gefunden, mit denen man MNL-Modelle schätzen kann aquí o aquí .

Ich habe einige andere Schriften über das "Rollen" Ihrer eigenen MLE-Funktion gefunden aquí . Wie auch immer, nach meinen Recherchen sind alle diese Funktionen und Pakete auf die interne optim Funktion.

In meinen Benchmark-Tests, optim ist der Engpass. Bei Verwendung eines simulierten Datensatzes mit ~16000 Beobachtungen und 7 Parametern benötigt R auf meinem Rechner etwa 90 Sekunden. Das entsprechende Modell in Biogeme dauert ~10 Sekunden. Ein Kollege, der seinen eigenen Code in Ox meldet etwa 4 Sekunden für dasselbe Modell.

Hat jemand Erfahrung mit dem Schreiben einer eigenen MLE-Funktion oder kann mich auf etwas hinweisen, das über die Standardfunktion hinaus optimiert ist? optim Funktion (kein Wortspiel beabsichtigt)?

Wenn jemand den R-Code zur Erstellung des Modells haben möchte, kann er ihn mir gerne zur Verfügung stellen. Ich habe ihn nicht zur Verfügung gestellt, da er nicht direkt für das Problem der Optimierung der optim Funktion und um Platz zu sparen...

EDIT: Vielen Dank an alle für ihre Gedanken. Ausgehend von einer Vielzahl von Kommentaren unten konnten wir R für kompliziertere Modelle in die gleiche Größenordnung wie Biogeme bringen, und R war sogar schneller für mehrere kleinere / einfachere Modelle, die wir ausgeführt haben. Ich denke, die langfristige Lösung für dieses Problem wird das Schreiben einer separaten Maximierungsfunktion beinhalten, die sich auf eine Fortran- oder C-Bibliothek stützt, aber ich bin sicherlich offen für andere Ansätze.

Answer 1

Schon mit der Funktion nlm() ausprobiert? Ich weiß nicht, ob es viel schneller ist, aber es verbessert die Geschwindigkeit. Überprüfen Sie auch die Optionen. optim verwendet als Standard einen langsamen Algorithmus. Sie können eine > 5-fache Beschleunigung erreichen, wenn Sie den Quasi-Newton-Algorithmus (method="BFGS") anstelle des Standardalgorithmus verwenden. Wenn Ihnen die letzten Ziffern nicht so wichtig sind, können Sie auch die Toleranzstufen von nlm() höher setzen, um zusätzliche Geschwindigkeit zu erreichen.

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

Answer 2

Ich bin der Autor des R-Pakets optimParallel die in Ihrem Fall hilfreich sein könnten. Das Paket bietet parallele Versionen der gradientenbasierten Optimierungsmethoden von optim() . Die Hauptfunktion des Pakets ist optimParallel() die dieselbe Verwendung und Ausgabe hat wie optim() . Verwendung von optimParallel() kann die Optimierungszeiten erheblich verkürzen, wie die folgende Abbildung zeigt ( p ist die Anzahl der Parameter). Siehe https://cran.r-project.org/package=optimParallel y http://arxiv.org/abs/1804.11058 für weitere Informationen.

Answer 3

Haben Sie das Material auf dem CRAN-Aufgabenansicht für die Optimierung ?

Answer 4

FWIW, ich habe das in C-ish gemacht, mit OPTIF9. Schneller geht's kaum. Es gibt viele Möglichkeiten, etwas langsamer zu machen, z. B. mit einem Interpreter wie R.

Hinzugefügt: Aus den Kommentaren geht hervor, dass OPTIF9 als Optimierungsmaschine verwendet wird. Das bedeutet, dass höchstwahrscheinlich der größte Teil der Zeit für die Auswertung der Zielfunktion in R verwendet wird. Es ist zwar möglich, dass C-Funktionen für einige der Operationen darunter verwendet werden, aber es gibt immer noch Interpreter-Overhead. Es gibt eine schnelle Methode, um herauszufinden, welche Codezeilen und Funktionsaufrufe in R den größten Teil der Zeit beanspruchen, und zwar indem man die Berechnung mit der Escape-Taste anhält und den Stapel untersucht. Wenn eine Anweisung X % der Zeit kostet, befindet sie sich zu X % der Zeit auf dem Stack. Sie werden feststellen, dass es Operationen gibt, die nicht nach C gehen, aber dort liegen sollten. Jeder Beschleunigungsfaktor, den Sie auf diese Weise erzielen, bleibt erhalten, wenn Sie einen Weg finden, die R-Ausführung zu parallelisieren.