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

[feedinglight]

Это среднеквадратичное?

[Andrey_Kireew]

нет, просто сумма квадратов

[feedinglight]

Настало время учиться адекватно оценивать погрешности..

И зачем надо было извращать алгоритм, который без изменений принимает на вход данный функционал? Может быть там ещё ошибки наляпаны, но теперь, хотя бы при $\alpha = 0$ и $\alpha=\infty$, вроде правильно (лень проверять в других точках).

Исходная задача наверно была не просто минимизировать функционал, а найти минимум, близкий к $\bold x = 0$, если сначала брать большое $\alpha$, а потом уменьшать. Не проще тогда было задать начальное приближение $\bold x = 0$?

[Andrey_Kireew]

feedinglight в сообщении #1394738 писал(а):

И зачем надо было извращать алгоритм, который без изменений принимает на вход данный функционал? Может быть там ещё ошибки наляпаны, но теперь ... вроде правильно (лень проверять в других точках).

Видимо, нелегко давать советы, и находить чужие ошибки, не имея элементарных представлений о предмете обсуждения. Проигнорировав вопрос

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

feedinglight в сообщении #1394627 писал(а):

та альфа по идее должна просто попасть в матрицу $\bold J$ и никуда не вылазить из неё,

это почему же?

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

[feedinglight]

Интересно было бы увидеть по-больше численных результатов

[vpb]

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

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

Andrey_Kireew
Вы сами весьма неважно объясняете, поэтому не стоит над слушателями иронизировать и т.д. У feedinglight здравый смысл вполне присутствует (хотя пока, кажется, знаний еще не много).

Объясним (для Евгений Машеров и feedinglight; впрочем, для первого большая часть нижеследующего тривиальна), какой именно вариант метода Гаусса-Ньютона используется. Обычная формулировка содержится в
Ортега, Рейнболдт, Итерационные методы решения систем со многими неизвестными, параграф 8.5; Дэннис, Шнабель, Численные методы безусловной оптимизации и решения нелинейных уравнений, пар. 10.2.

Рассмотрим отображение $F\colon{\mathbb R}^m\longrightarrow {\mathbb R}^n$, класса $C^1$. (ТС называет $\delta=(F_1(x),\ldots,F_n(x))^t$ "вектором невязок", соответствующим данному $x$).

Мы хотим найти $x\in{\mathbb R}^m$ такое, что $F(x)=\overline0$, или, по крайней мере, такое, на котором функция $r(x)=\sum_{i=1}^nF_i^2(x)$ ($\xi$, в обозначениях ТС) достигает локального минимума.

Пусть $F'(x)$ --- производная отображения $F$ в точке $x$, т.е.

$$ F'(x)=J(x)=\left( \frac{\partial F_i}{\partial x_j} (x)\right) $$

--- матрица Якоби (размера $n\times m$).
Обычный метод Ньютона есть

$$ x^{k+1}=x^k-J(x^k)^{-1} F(x^k) $$

(для случая, когда $m=n$, а матрица $J(x)$ невырождена, в том числе в точке $x_\ast$, являющейся искомым решением уравнения $F(x_\ast)=0$). (Замечание. Метод Ньютона для системы $F(x)=0$, и метод Ньютона для минимизации функции $r(x)$ --- это несколько разные вещи).

Итерационный процесс

$$ x^{k+1}=x^k-\lambda_k J(x^k)^{-1} F(x^k) $$

--- демпфированный метод Ньютона. $\lambda_k$ выбираются так, чтобы обеспечить убывание $r(x^{k+1})<r(x^k)$ (например).

-- 26.05.2019, 19:56 --

В общем же случае (когда $m\ne n$ или $J(x)$ вырождено) формула модифицируется как

$$ x^{k+1}=x^k-\lambda_k J(x^k)^+\ F(x^k), $$

где $A^+$ --- псевдообратная по Муру-Пенроузу матрицы $A$. Это (демпфированный) метод Гаусса-Ньютона.
Если $A$ --- $n\times m$ матрица, причем $m\leq n$ и $\operatorname{rk}(A)=m$, то $A^+=(A^tA)^{-1}A^t$. В этом случае формула приобретает вид

$$ x^{k+1}=x^k-\lambda_k (J(x^k)^tJ(x^k))^{-1}J(x^k)^t\ F(x^k).$$

Вообще говоря, взятие псевдообратной --- неустойчивая операция. Причина в том, что $A^tA$ может быть вырождена. Поэтому прибегают к такой регуляризации:

$$ x^{k+1}=x^k-\lambda_k (J(x^k)^tJ(x^k)+\alpha_kE_m)^{-1}J(x^k)^t\ F(x^k), $$

где $\alpha_k>0$. Это называется "метод Левенберга-Марквардта". Формула из стартового поста ТС --- это как раз формула метода ЛМ. (Замечание. Разные названия, связанные с методом Ньютона, часто употребляются одни вместо других. Тут есть путаница).

ТС, как я понял, использует следующую регуляризацию ГН (в исправленном посте). Возьмем какое-либо $\beta>0$. Рассмотрим такое отображение

$$ \widetilde F=\widetilde F_\beta\colon{\mathbb R}^m\longrightarrow {\mathbb R}^{n+m}, \qquad \widetilde F_\beta(x)=(F(x), \beta x)^t\,.$$

Также рассмотрим

$$\widetilde r_\beta(x)=|\widetilde F_\beta(x)|^2= \sum_{i=1}^n F_i^2(x) +\beta^2|x|^2.$$

Ясно, что когда $\beta\longrightarrow0$, то точка минимума для $\widetilde r_\beta(x)$ стремится к таковой для $r(x)$. С другой стороны, соответствующая матрица Якоби есть

$$ \widetilde J_\beta(x)=\begin{pmatrix} J(x) \\ \beta E_m \end{pmatrix} $$

и тем самым имеет ранг $m$. Поэтому можно рассмотреть демпфированную итерацию Гаусса-Ньютона для $\widetilde F_\beta(x)$. Она имеет вид

$$ x^{k+1}=x^k-\lambda_k (\widetilde J_\beta(x^k)^t\widetilde J_\beta(x^k))^{-1}\widetilde J_\beta(x^k)^t\ \widetilde F_\beta(x^k).$$

После небольших преобразований последняя формула переписывается как

$$ x^{k+1}=x^k-\lambda_k (J(x^k)^tJ(x^k)+\beta^2E_m)^{-1}(J(x^k)^t F(x^k)+\beta^2x).$$

Обозначая $\beta^2=\alpha$, получим ту же формулу, что и ТС.

[Andrey_Kireew]

В дополнение к обзору vpb приведу полноценный алгоритм Ньютона для МНК:

$$ \bold x_{k+1}=\bold x_k-\lambda_k (\bold J_k^T \bold J_k+\sum\limits_{i=1}^{m}\delta_i \bold H_i+\alpha \bold E_m)^{-1}(\bold J_k^T \bold \delta +\alpha \bold [x-x_0]). \qquad (1) $$

Нетрудно заметить, что все описанные ранее алгоритмы являются частным случаем (1).

P.S.: алгоритм Левенберга-Марквардта работает немного по другому: у него нет настраиваемого $\lambda_k$, вместо этого используется настраиваемый $\alpha$ (в методе Ньютона-Гауса он фиксирован).
Вообще, это эвристический алгоритм, не имеющий строгого обоснования. Обычно приводится такое объяснение: при большом $\alpha $ алгоритм приближается к градиентному поиску, становится устойчивее, но сходится медленнее. Иногда он представляется в такой модификации

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

Я такой как то пробовал. Кажется работает лучше, чем с единичной диагональю. Но в таких делах всё зависит от конкретной задачи.
Ну, и самый очевидный недостаток - подбор $\alpha $ требует многократного обращения матрицы на каждой итерации, что весьма накладно.

[vpb]

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

полноценный

Поясним для посторонего читателя, что в данном случае слово "полноценный" /"неполноценный" не несет эмоциональной окраски, а означает класс алгоритмов с полным или неполным гессианом (с неполным --- это Гаусс-Ньютон, Левенберг-Марквардт.). Неполноценные в Деннис-Шнабеле трактуются в пар. 10.2, а полноценные в 10.3. Те и другие имеют свои преимущества и недостатки.

Дисклеймер: сам я по оптимизации не являюсь особым специалистом, если что.

[Andrey_Kireew]

Ответ на мой вопрос был найден в книге, любезно предоставленной vpb. Как ни странно, но метод Ньютона-Гаусса в общем случае не обладает свойством локальной сходимости. Об этом я и не подозревал. Для этого необходимо выполнение ряда специальных условий. В общем, с приближением к оптимальному решению сходимость у него замедляется. При некоторой остаточной невязке оптимизация может вообще остановится. Собственно это у меня и наблюдается. До предела вычислительной точности ещё очень далеко, а итеративный процесс уже перешел в установившийся режим, невязка практически не уменьшается.
Выясняется, что дело тут совсем не в регуляризации. Просто в конкретной задаче введение дополнительного слагаемого в градиент несколько изменило форму целевой функции, что и обеспечило условие локальной сходимости. На других примерах было обнаружено, что в общем случае это не решение проблемы. Даже с изменённым градиентом алгоритм так же часто останавливается, не достигнув оптимума.

[feedinglight]

vpb
Спасибо за объяснение. Да, случай $\alpha = 0$ я рассматривал не совсем корректно. Нужно чтобы матрица обращалась.

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

После небольших преобразований последняя формула переписывается как

$$ x^{k+1}=x^k-\lambda_k (J(x^k)^tJ(x^k)+\beta^2E_m)^{-1}(J(x^k)^t F(x^k)+\beta^2x).$$

Ну да, вроде правильно. То есть итоговая схема ТС - это просто схема Ньютона Гаусса, применённая к заданному функционалу $\xi=\bold \delta^T \bold \delta+\alpha \bold x^T\bold x.$. Значит то, что я назвал "извращением" - оптимизированный для реализации вариант записи.

Andrey_Kireew
Множители Лагранжа не пробовали?

[vpb]

feedinglight в сообщении #1395676 писал(а):

Множители Лагранжа не пробовали?

Э... непонятно, какое отношение к вопросу имеют множители Лагранжа. В теме же речь о безусловной оптимизации, а множители Лагранжа используются в условной. Или есть какие-то другие множители Лагранжа, о которых я не в курсе ?

[Andrey_Kireew]

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

В теме же речь о безусловной оптимизации, а множители Лагранжа используются в условной

просто эта тема является логическим продолжением предыдущей, в обсуждении которой feedinglight тоже участвовал. Поэтому он знает, что тут действительно оптимизация условная. Там есть строгое ограничение на положительные коэффициенты. Хорошее решение предложил Евгений Машеров - сделать перепараметризацию $y=log(x)$.
От множителей Лагранжа я пока особой пользы не вижу. Но я с ними и не знаком. Конечно, может в этом что то и есть.

[vpb]

Понятно. Я думал, что эта тема отдельная.

[Andrey_Kireew]

Всё таки удалось реализовать полноценный алгоритм Ньютона. Быстро выяснилось, что причина "застревания" алгоритма вовсе не в игнорировании слагаемого со вторыми производными в матрице Гёссе. После ряда экспериментов, для некоторых начальных приближений, "застревание" обнаружилось и для полноценного алгоритма. Сначала было предположение, что алгоритм попадает в локальные минимумы. Но в последствие выяснилось, что это не так. В действительности, причина в нарушении условия сходимости алгоритма Ньютона (его можно найти в книге, любезно предоставленной vpb, на стр. 141):

$$\bold \nabla^T\bold H^{-1}\bold \nabla>0,\qquad (1)$$

где $\bold \nabla =\bold J^T \bold \delta +\alpha \bold x,$ $\bold H=\bold J^T \bold J+\alpha \bold E. $
Анализ изменения критерия (1) от итерации к итерации показал, что "застревание" происходит происходит сразу же после его нарушения.

Решение проблемы у меня пока такое: на каждой итерации выполняется проверка условия (1). Если оно выполняется, то осуществляется шаг в Ньютоновском направлении, с измельчением его величины до тех пор, пока не уменьшится невязка. Если же (1) не выполняется, то выполняется шаг в направлении антиградиента. Кажется, пока работает. Ещё есть идея - увеличивать $\alpha$, до тех пор, пока (1) не будет выполнено. Наверное, второй вариант будет лучше первого.

Но тут есть проблема. Условие (1) можно проверить только в полноценном алгоритме Ньютона, а он очень медленный. В методе Гаусса-Ньютона этот способ уже не помогает "застревания" иногда происходят. Можно, конечно модифицировать (1) для этого случая так

$$\bold \nabla^T\bold H^{-1}\bold \nabla>C,\qquad (2)$$

но тогда возникают проблемы с выбором константы C.

[vpb]

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

книге, любезно предоставленной vpb, на стр. 141):

Пожалуйста, на здоровье. (Для третьих лиц: имеется в виду Дэннис-Шнабель).

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

Гёссе

Правильно Гессе, без ё.

Непонятно, как может нарушаться условие (1). Матрица $H=J^TJ+\alpha E$ положительно определена (и симметрична), поэтому и $H^{-1}$ положительно определена, для любой матрицы $J$ и любого $\alpha>0$, откуда $v^T H^{-1} v>0$ для любого ненулевого вектора $v$.

-- 01.06.2019, 02:55 --

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