Остаточная ошибка метода Гаусса-Ньютона

[Andrey_Kireew]

vpb в сообщении #1396985 писал(а):

Но $H'$, и $H'+\alpha E$ тоже, где $H'$ обозначает полный гессиан, может не быть положительно определенной

Возможно, действительно, всё дело в добавочном члене. Сейчас проверил: минимальные собственные числа полной $ \bold H $ имеют порядок $10^{-5}$, а без добавочного члена - уже $10^{-20}$. Понятно, что речь идет об отрицательных числах.

[vpb]

У $J^TJ+\alpha E$ при $\alpha\geq0$ отрицательных собственных значений вообще быть не может, разве что из-за ошибок округления.

[Andrey_Kireew]

vpb в сообщении #1397126 писал(а):

У $J^TJ+\alpha E$ при $\alpha\geq0$ отрицательных собственных значений вообще быть не может, разве что из-за ошибок округления

Да, похоже на ошибки округления. Тем не менее, исходный вопрос так и остаётся открытым.
Эксперименты показывают, что существует оптимальное значение $\alpha$, при котором достигается минимальная невязка. Причём, она далека от предела, определяемого ошибками округления. Кривая сходимости выглядит так:

Создаётся впечатление, что алгоритм во что то "упирается".

Но опытным путём, был найден такой способ оптимизации: сначала оптимизация выполняется с большим $\alpha$, до тех пор, пока процесс не установится. Затем $\alpha$ уменьшается. При этом скорость сходимости (непонятно почему) возрастает. После того, как уменьшение невязки существенно замедлится, $\alpha$ снова уменьшается, и т.д.
Таким образом удаётся довести невязку до предела вычислительной погрешности.
Вот пример графика сходимости, для той же модели и того же начального приближения, что и выше. Параметр $\alpha$ уменьшался в 10 раз через каждые 100 итераций:

Значение $\sim 10^{-32}$ - это как раз погрешность вычисления невязки.

Если же сразу поставить маленький $\alpha$, алгоритм "упрётся" в $\sim 10^{-7}$ и на этом остановится.

[vpb]

Andrey_Kireew
Толком я вам, наверное, ничего объяснить не смогу. Т.е. у меня нет теоретического объяснения, почему Гаусс-Ньютон так себя ведет. Собственно, вообще он не входит в мои интересы.

Козьма Прутков писал(а):

Многие вещи непонятны нам не потому, чтобы наши понятия слабы, но потому, что сии вещи не входят в круг наших понятий.

Я, видите ли, в силу специфики своих задач предпочитаю Ньютона с полным гессианом.

Пословица писал(а):

Девяносто пять процентов мужчин любят полных женщин. И лишь пять процентов --- очень полных.

По поводу полного Ньютона: такое поведение, по моим соображениям, тоже вполне возможно. Надо, однако, уточнить.

Есть, по существу, два варианта метода Ньютона: с линейным поиском, и с доверительной областью. Линейный поиск, по мне, это ненадежно. "Чистое" ньютоновское направление может чёрти куда указывать, и допустимый шаг вдоль него может быть очень маленьким, что создает впечатление, что застряли. Можно добавить регуляризующее слагаемое $\alpha E$, получается алгоритм

$$ x_{n+1}=x_n-\mu_n(Hf(x_n)+\alpha_n E)^{-1} \nabla f(x_n) $$

где $\alpha_n$ как-то (непонятно, как; в разных случаях по разному) выбирается, а $\mu_n$ берется из условия убывания (Армихо и др.). (В этой записи подразумевается, что мы рассматриваем задачу минимизации функции $f$, а $\nabla f$ и $Hf$ --- ее градиент и гессиан).
(продолжение следует).

-- 03.06.2019, 01:06 --

В Дэннис-Шнабеле, емнип, $\alpha_n$ выбирается так, чтобы $Hf(x_n)+\alpha_n E$ было положительно определено и достаточно хорошо обусловлено, скажем число обусловленности было не более чем порядка $\operatorname{macheps}^{-1}$. Ну а $\mu_n$ подбирается потом.

Точно я сейчас не могу сформулировать, но, в общем, при таком подходе как раз возможно описанное вами поведение алгоритма. Особенно в случае, когда гессиан оказывается плохо обусловлен, а среди его собственных значений есть и маленькие, и большие, и промежуточные ($10^{-5}$, $1$, $10^5$). Размышлять, как именно это может случиться, я, извините, не хочу, потому что думаю так: зачем размышлять об особенностях плохого алгоритма, когда есть хороший, надежный во всех случаях ? Под надежным я понимаю второй вариант метода Ньютона, а именно с доверительной областью.

(продолжение следует)

-- 03.06.2019, 01:18 --

В общем, того же и Вам рекомендую. А прочитать можно (а может, Вы всё это и сами знаете...) можно: Дэннис-Шнабель, глава 6; Nocedal, Wright, Numerical optimization, и целая книжка Conn, Gould, Toint, Trust region methods.

(кстати: сведения о первых двух книжках я почерпнул именно с нашего форума. Впоследствии, читая литературу, узнал, что это классика).
(ну вот, вроде, всё)

[Andrey_Kireew]

Решил внять Вашему vpb совету, и опробовать алгоритм с регулируемым $\alpha$. Вот такой:

$$ \bold x_{k+1}=\bold x_k-(\bold J_k^T \bold J_k+\alpha \bold E_m)^{-1}\bold J_k^T \bold \delta. \qquad (1) $$

Стратегию изменения $\alpha$, для начала, выбрал самую простую, изначально предложенную Марквардтом [Демиденко Е.З. Линейная и нелинейная регрессии. М.: Финансы и статистика, 1981. -С. 254.], а именно - некоторое начальное $\alpha$ на каждой итерации уменьшается в несколько раз. Если, после этого функционал не убывает - $\alpha$ увеличивается, до тех пор, пока не будет достигнуто убывание. В указанном выше источнике коэффициент уменьшения и коэффициент увеличения $\alpha$ равны 10, а начальное значение берётся 0.01. Вычисления показывают, что при этих значениях, на моей задаче, наблюдается такое же "застревание" как и в алгоритме с дроблением шага. При этом $\alpha$ стремится к бесконечности и алгоритм останавливается, не достигнув минимума целевой функции.
Я использовал коэффициент дробления $1,5$ и коэффициент увеличения $2$, начальное значение оставил таким же (оказалось, что оно мало на что влияет). В алгоритма Марквардта используется пропорциональный регуляризатор, т.е. вместо единичной матрицы, в алгоритме используется диагональ $\bold J_k^T \bold J_k$. Встречал в разных источниках различные описания преимуществ данного подхода. Однако, у меня, всё получилось, с точностью до наоборот. Пропорциональный регуляризатор оказался хуже. Алгоритм с ним сходится быстрее только на начальном этапе, а потом, скорость сходимости сильно замедляется. Хуже всего - с ним алгоритм более склонен "застреванию". Его работа так же существенно зависит от выбора начального $\alpha$.
С прямым же регуляризатором удалось добиться более-менее устойчивой работы алгоритма. Вычисления показывают, что теперь, он сходится от самых разных начальных приближений. Кривая сходимости всегда примерно одинаковая, и имеет следующий вид


При достижении предела вычислительной точности $\alpha$ очень быстро возрастает до $\infty$. Поэтому, мне показалось очень удобным использовать следующий критерий останова:

$$\alpha>10^{100}$$

.
Так же выяснился один нюанс - точность процедуры обращения матрицы здесь очень важна. У меня используется весьма трудоёмкая процедура псевдообращения, на основе сингулярного разложения. Попытка замены этой процедуры прямым алгоритмом решения СЛАУ, на основе LU разложения, приводит к полной потере работоспособности. Наверное, это актуально только для больших, плохо обусловленных матриц.

Разумеется, надо будет попробовать более совершенные стратегии подбора $\alpha$, описанные в книге Дэннис-Шнабель, а так же - использовать полную матрицу Гессе.
Вообще, конечно, сходится полученный алгоритм не так быстро, как хотелось бы. Тем более, что задача была выбрана специально небольшая.
Но это уже кое что! Надеюсь, в дальнейшем получиться его улучшить.
Большое спасибо Вам vpb, за помощь

[Andrey_Kireew]

Кажется, причина интересного явления

Andrey_Kireew в сообщении #1397342 писал(а):

[quote="vpb в [url=http://dxdy.ru/post1397126.html#p1397126] Кривая сходимости выглядит так:

Создаётся впечатление, что алгоритм во что то "упирается".

обнаружена.

Оказалось, что такое же поведение, иногда, наблюдается и у градиентного алгоритма. Проблема здесь в ошибке вычисления градиента.
При приближении к оптимальному решению градиент становится мал и ошибка его вычисления - существенной. Обозначенная проблема возникает, когда направление вычисленного градиента, из за ошибок, перестаёт быть направлением спуска. Целевая функция в этом направлении только возрастает. Все попытки подобрать оптимальный шаг (путём измельчения или другими способами) приводят к исчезающе малым шагам. На новой итерации, из-за очень малого шага, градиент вычисляется почти в той же самой точке, и ситуация повторяется. Так как машинная точность ограничена, шаг быстро становится равным нулю и возникает это "застревание".

[vpb]

Andrey_Kireew
Не знаю что сказать ... Посоветую только овладеть методом доверительной области. Только у меня самого недавно с этим методом случился афронт, или как там оно называется, в общем, неожиданная неудачка. Думал я долго (месяца 4), что я понимаю, что это такое, но вдруг обнаружил в своих понятиях об этом предмете большой пробел, и пока не разобрался.