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

[Andrey_Kireew]

Имеется следующий алгоритм многомерной оптимизации:

$$\bold x_{n+1}=\bold x_{n+1}-\mu (\bold J^T\bold J^T+\alpha \bold E)^{-1} \bold J^T \bold \delta, $$

где $\bold J$ - матрица Якоби, $\bold \delta$ - вектор невязок, $\mu \in (0, 1] $ - меняющийся шаг, $\alpha $ - параметр регуляризации, $\bold E$ - единичная диагональная матрица.
Оптимизируемый функционал имеет вид:

$$\xi=\bold \delta^T \bold \delta+\alpha \bold x^T\bold x.$$

С уменьшением $\alpha$, как и следовало ожидать, остаточная невязка $\xi_{min}$ уменьшается. Но это наблюдается лишь до некоторого $\alpha_{min}$. При дальнейшем уменьшении $\alpha$ остаточная невязка увеличивается, и при $\alpha=0$ она существенно больше, чем при $\alpha_{min}$.

В общем, до конца алгоритм сходиться не хочет. До eps ещё очень далеко. В чём может быть дело?

[feedinglight]

При $\alpha = 0$ получится метод Ньютона, в котором вместо $\bold J$ почему-то взята $\bold J^T$. Что это за метод такой?

[Евгений Машеров]

Там точно как приведено? Транспонированная на транспонированную?

[Andrey_Kireew]

Понятно, что опечатка. Правильно так:

$$\bold x_{n+1}=\bold x_{n+1}-\mu (\bold J^T\bold J+\alpha \bold E)^{-1} \bold J^T \bold \delta. $$

Строго говоря, методом Ньютона это является только в случае нулевой остаточной невязки. Т.е. при точном решении СНАУ. Если невязка не нулевая, например, когда $\bold J$ прямоугольная (а у меня как раз так и есть), то направление оптимизации отличается от ньютоновского, так как не учитываются вторые производные. Поэтому он и называется метод Ньютона-Гаусса.

[Евгений Машеров]

В качестве простейшего объяснения - матрица $\bold J$ неполного ранга или близка к матрице неполного ранга, производные, они часто оказываются сильно скоррелированными (просто по формуле производной сложной функции - общие сомножители появляются), "мультиколлинеарность", как это зовут применительно к регрессии. Матрица $\bold J^T\bold J^T+\alpha \bold E $ при обнулении альфы получается вырожденной или близкой к вырожденной. В силу конечной точности вычислений при попытке ея обращения

(Оффтоп)

в истинную веру

вместо честного деления на ноль и "авоста" получаем игру ошибок округления (и рост элементов обращённой до бесконечности, но при умножении на $\bold J^T$ большие сомножители гасятся маленькими, и результат выглядит правдоподобным, но никакого отношения к оптимизации уже не имеет, это именно шуточки округления. Прибавка единичной матрицы это предотвращает. Здесь хорошо проясняет сделать сингулярное разложение для $\bold J$, сократить и увидеть, что происходит с сингулярными числами в решении и как работает

(Оффтоп)

спецподразделение

альфа.

[Andrey_Kireew]

Насколько я понимаю, обусловленность $\bold J $ влияет только на ошибки $\bold x$, которые определяются так:

$$|\delta \bold x|=cond(\bold J)\cdot\xi,$$

но сама невязка $\xi$ не должна зависеть от $cond(\bold J)$. Она всё равно должна уменьшаться до некоторого минимального значения, а этого не происходит. Я специально использовал точную модель, для которой $\xi\approx 0$ (порядка $10^{-34}$), но невязка уменьшается только до $10^{-8}$-$10^{-10}$. Непонятно, почему она не уменьшается дальше?

[Евгений Машеров]

Величина, достигнутая при оптимизации, вообще от альфы не зависит. Она влияет только на скорость движения к оптимуму. При очень большой альфе получается движение по градиенту мелкими шагами. То есть медленнее, но надёжнее.
Двойная точность это погрешность порядка $10^{-15}$, по числу разрядов мантиссы, а если брать не единичную операцию, а сложное вычисление - $10^{-10}$, откуда $10^{-34}$?

[Andrey_Kireew]

это абсолютная невязка, она может быть и $10^{-300}$

[Евгений Машеров]

Без знания вида задачи не могу говорить более определённо, но у меня впечатление, что просто вышли на максимум точности. Во всяком случае, значение альфы тут не столь важно. Нулевая альфа - пытаемся одним прыжком напрыгнуть на оптимум, но есть шанс из-за численной неустойчивости быть заброшенными далеко от него. Большая альфа - получается почти в точности градиентный поиск помаленьку. А оптимум тот же самый. Просто быстрее или медленнее приближаемся.

[Andrey_Kireew]

Евгений Машеров в сообщении #1394619 писал(а):

но у меня впечатление, что просто вышли на максимум точности

В этом и весь вопрос. Если взять начальное приближение, равное точному решению, то невязка получается $10^{-34}$. Если взять очень близко к истинному решению - получается уже $10^{-34}$. Если начальное приближение достаточно далеко, то получается примерно $10^{-8}$. Создаётся впечатление, что при достижении некоторой невязки, например $10^{-8}$, процесс оптимизации останавливается, ну или очень сильно замедляется.

[feedinglight]

Если минимизируется $\xi=\bold \delta^T \bold \delta+(\bold {\sqrt{\alpha} x})^T(\bold {\sqrt{\alpha} x})$, то почему прибавляется $\alpha \bold E$. Эта альфа по идее должна просто попасть в матрицу $\bold J$ и никуда не вылазить из неё, или это другая альфа? При $\alpha=\infty$ должно получаться $\bold x = 0$. Не понятно.

-- 22.05.2019, 21:07 --

а невязка, если это сумма квадратов, совсем плохая (

[Andrey_Kireew]

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

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

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

[Евгений Машеров]

И можно прояснить

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

Оптимизируемый функционал имеет вид:

$$\xi=\bold \delta^T \bold \delta+\alpha \bold x^T\bold x.$$

?

[Andrey_Kireew]

что конкретно прояснить?

[Andrey_Kireew]

Благодаря намёку Евгений Машеров я нашел ошибку в алгоритме. Правильно должно быть так:

$$\bold x_{n+1}=\bold x_{n}-\mu (\bold J^T\bold J+\alpha \bold E)^{-1} (\bold J^T \bold \delta+\alpha \bold x_{n}). $$

дело, видимо, было в последнем слагаемом.

Теперь, кажется, всё сходится до предела ошибок вычислений: