Как обрабатывать граничные ограничения при использовании `nls.lm` в R

Я спросил этот вопрос некоторое время назад. Я не уверен, должен ли я опубликовать это как ответ или новый вопрос. У меня нет ответа, но я "решил" проблему, применив алгоритм Левенберга-Марквардта, используя nls.lm в R, и когда решение находится на границе, я запускаю алгоритм рефлексии в области доверия (TRR, реализованный в R), чтобы отойти от него. Теперь у меня есть новые вопросы.

Из моего опыта, таким образом, программа достигает оптимальной и не очень чувствительна к стартовым значениям. Но это всего лишь практический метод, чтобы отвлечься от проблем, с которыми я столкнулся, используя nls.lm, а также другие функции оптимизации в R. Я хотел бы знать, почему nls.lm ведет себя таким образом для задач оптимизации с граничными ограничениями и как обрабатывать граничных ограничений при использовании nls.lm на практике.

Далее я привел пример, иллюстрирующий два вопроса, используя nls.lm.

Он чувствителен к стартовым значениям.
Он останавливается, когда какой-то параметр достигает границы.

Воспроизводимый пример: Фокусный набор D

library(devtools)
install_github("KineticEval","zhenglei-gao")
library(KineticEval)
data(FOCUS2006D)
km <- mkinmod.full(parent=list(type="SFO",M0 = list(ini   = 0.1,fixed = 0,lower = 0.0,upper =Inf),to="m1"),m1=list(type="SFO"),data=FOCUS2006D)
system.time(Fit.TRR <- KinEval(km,evalMethod = 'NLLS',optimMethod = 'TRR'))
system.time(Fit.LM <- KinEval(km,evalMethod = 'NLLS',optimMethod = 'LM',ctr=kingui.control(runTRR=FALSE)))
compare_multi_kinmod(km,rbind(Fit.TRR$par,Fit.LM$par))
dev.print(jpeg,"LMvsTRR.jpeg",width=480)

Дифференциальные уравнения, описывающие модель/систему:

"d_parent = - k_parent * parent"                                                 
"d_m1 = - k_m1 * m1 + k_parent * f_parent_to_m1 * parent"

В графе слева находится модель с начальными значениями, а в середине - модель с использованием "TRR" (аналогичная алгоритму в функции Matlab lsqnonlin), справа - установленная модель с использованием "LM" с nls.lm. Изучив установленные параметры (Fit.LM$par), вы обнаружите, что один установленный параметр (f_parent_to_m1) находится на границе 1. Если я изменил начальное значение для одного параметра M0_parent от 0,1 до 100, то я получил те же результаты, используя nls.lm и lsqnonlin. У меня много таких случаев.

newpars <- rbind(Fit.TRR$par,Fit.LM$par)
rownames(newpars)<- c("TRR(lsqnonlin)","LM(nls.lm)")
newpars

                   M0_parent   k_parent        k_m1 f_parent_to_m1
TRR(lsqnonlin)  99.59848 0.09869773 0.005260654       0.514476
LM(nls.lm)      84.79150 0.06352110 0.014783294       1.000000

За исключением вышеперечисленных проблем, часто бывает, что гессиан, возвращаемый nls.lm, не является инвертируемым (особенно когда некоторые параметры находятся на границе), поэтому я не могу получить оценку ковариационной матрицы. С другой стороны, алгоритм "TRR" (в Matlab) почти всегда дает оценку, вычисляя якобиан в точке решения. Я думаю, что это полезно, но я также уверен, что алгоритмы оптимизации R (те, которые я пробовал) не сделали этого по какой-то причине. Я хотел бы знать, ошибаюсь ли я, используя метод Matlab для вычисления матрицы ковариации, чтобы получить стандартную ошибку для оценок параметров.

Последнее замечание, которое я утверждал в предыдущем сообщении

Ответ 1

Во-первых, я не эксперт по Matlab и Optimization и никогда не использовал R.

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

LM слегка усилен подход Гауза-Ньютона - для задач с несколькими локальными минимумами он очень чувствителен к начальным состояниям. Включение границ обычно генерирует больше этих минимумов.

TRR сродни LM, но более устойчив. Он имеет лучшие возможности для "выпрыгивания" из плохих локальных минимумов. Вполне возможно, что он будет вести себя лучше, но хуже, чем LM. На самом деле объяснить, почему очень сложно. Вам нужно будет подробно изучить алгоритмы и посмотреть, как они себя ведут в этой ситуации.

Я не могу объяснить разницу между реализацией Matlab и R, но есть несколько расширений TRR, которые, возможно, используются Matlab, а R - нет. Ваш подход к использованию LM и TRR чередуется сходится лучше, чем TRR?

Ответ 2

Используя пакет mkin, вы можете найти параметры, используя алгоритм "Порт" (который также является своего рода алгоритмом TRR, насколько я могу судить по его документации) или алгоритм "Marq", который использует nls.lm в фоновом режиме. Затем вы можете использовать "обычные" начальные значения или "плохие" начальные значения.

library(mkin)
packageVersion("mkin")

Последняя версия mkin может значительно ускорить процесс, поскольку они компилируют модели из автоматически созданного кода C, если компилятор доступен в вашей системе (например, у вас установлен r-base-dev на Debian/Ubuntu или Rtools в Windows).

Это определяет модель:

m <- mkinmod(parent = mkinsub("SFO", "m1"),
             m1 = mkinsub("SFO"),
             use_of_ff = "max")

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

cat(m$diffs, sep = "\n")

Затем мы вписываемся в четыре варианта: Port и LM, с M0 или без него, фиксированным до 0,1:

f.Port = mkinfit(m, FOCUS_2006_D)
f.Port.M0 = mkinfit(m, FOCUS_2006_D, state.ini = c(parent = 0.1, m1 = 0))
f.LM = mkinfit(m, FOCUS_2006_D, method.modFit = "Marq")
f.LM.M0 = mkinfit(m, FOCUS_2006_D, state.ini = c(parent = 0.1, m1 = 0),
               method.modFit = "Marq")

Затем мы рассмотрим результаты:

results <- sapply(list(Port = f.Port, Port.M0 = f.Port.M0, LM = f.LM, LM.M0 = f.LM.M0),
  function(x) round(summary(x)$bpar[, "Estimate"], 5))

которые

                   Port  Port.M0       LM    LM.M0
parent_0       99.59848 99.59848 99.59848 39.52278
k_parent        0.09870  0.09870  0.09870  0.00000
k_m1            0.00526  0.00526  0.00526  0.00000
f_parent_to_m1  0.51448  0.51448  0.51448  1.00000

Итак, мы видим, что алгоритм Port находит лучшее решение (насколько мне известно) даже при плохих стартовых значениях. Проблема скорости, которая может иметься с более сложными моделями, облегчается с помощью автоматической генерации кода C.