Алгоритм Гаусса — Ньютона

Алгоритм Гаусса — Ньютона используется для решения задач нелинейным методом наименьших квадратов. Алгоритм является модификацией метода Ньютона для нахождения минимума функции. В отличие от метода Ньютона, алгоритм Гаусса — Ньютона может быть использован только для минимизации суммы квадратов, но его преимущество в том, что метод не требует вычисления вторых производных, что может оказаться существенной трудностью.

Приближение кривой со случайным шумом асимметричной моделью пика используя алгоритм Гаусса — Ньютона с переменным коэффициентом затухания α.
Сверху: необработанные данные и модель.
Снизу: ход нормализованной суммы квадратов ошибок.

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

Метод назван именами математиков Карла Фридриха Гаусса и Исаака Ньютона.

Описание

Если задано m функций r = (r₁, …, r_m) (часто называемых невязками) от n переменных β = (β₁, …, β_n), при m ≥ n. Алгоритм Гаусса — Ньютона итеративно находит значения переменных, которые минимизируют сумму квадратов

$S({\boldsymbol {\beta }})=\sum _{i=1}^{m}r_{i}^{2}({\boldsymbol {\beta }}).$

Начав с некоторого начального приближения ${\boldsymbol {\beta }}^{(0)}$, метод осуществляет итерации

${\boldsymbol {\beta }}^{(s+1)}={\boldsymbol {\beta }}^{(s)}-\left(\mathbf {J_{r}} ^{\mathsf {T}}\mathbf {J_{r}} \right)^{-1}\mathbf {J_{r}} ^{\mathsf {T}}\mathbf {r} ({\boldsymbol {\beta }}^{(s)})$

Здесь, если рассматривать r и β как вектор-столбцы, элементы матрицы Якоби равны

$(\mathbf {J_{r}} )_{ij}={\frac {\partial r_{i}({\boldsymbol {\beta }}^{(s)})}{\partial \beta _{j}}}$

а символ $^{\mathsf {T}}$ означает транспонирование матрицы.

Если m = n, итерации упрощаются до

${\boldsymbol {\beta }}^{(s+1)}={\boldsymbol {\beta }}^{(s)}-\left(\mathbf {J_{r}} \right)^{-1}\mathbf {r} ({\boldsymbol {\beta }}^{(s)})$,

что является прямым обобщением одномерного метода Ньютона.

При аппроксимации данных, где целью является поиск параметров β, таких, что заданная модель функций y = f(x, β) наилучшим образом приближает точки данных (x_i, y_i), функции r_i являются остаточными ошибками

$r_{i}({\boldsymbol {\beta }})=y_{i}-f(x_{i},{\boldsymbol {\beta }}).$

Тогда метод Гаусса — Ньютона можно выразить в терминах якобиана J_f функции f

${\boldsymbol {\beta }}^{(s+1)}={\boldsymbol {\beta }}^{(s)}+\left(\mathbf {J_{f}} ^{\mathsf {T}}\mathbf {J_{f}} \right)^{-1}\mathbf {J_{f}} ^{\mathsf {T}}\mathbf {r} ({\boldsymbol {\beta }}^{(s)}).$

Заметим, что $\left(\mathbf {J_{f}} ^{\mathsf {T}}\mathbf {J_{f}} \right)^{-1}\mathbf {J_{f}} ^{\mathsf {T}}$ является псевдообратной матрицей к $\mathbf {J_{f}}$.

Замечания

Требование m ≥ n в алгоритме необходимо, поскольку в другом случае матрица J_r^TJ_r не имеет обратной и нормальные уравнения нельзя решить (по меньшей мере однозначно).

Алгоритм Гаусса — Ньютона можно получить с помощью линейного приближения вектора функций r_i. Используя теорему Тейлора, мы можем для каждой итерации записать:

$\mathbf {r} ({\boldsymbol {\beta }})\approx \mathbf {r} ({\boldsymbol {\beta }}^{s})+\mathbf {J_{r}} ({\boldsymbol {\beta }}^{s})\Delta$,

где $\Delta ={\boldsymbol {\beta }}-{\boldsymbol {\beta }}^{s}$. Задача нахождения Δ, минимизирующего сумму квадратов в правой части, т.е.

$\mathbf {min} \|\mathbf {r} ({\boldsymbol {\beta }}^{s})+\mathbf {J_{r}} ({\boldsymbol {\beta }}^{s})\Delta \|_{2}^{2}$,

является линейной задачей нахождения наименьших квадратов, которую можно решить явно, что даёт нормальные уравнения.

Нормальные уравнения — это m линейных уравнений по неизвестным приращениям Δ. Уравнения могут быть решены за один шаг, если использовать разложение Холецкого, или, лучше, QR-разложение матрицы J_r. Для больших систем итеративный метод может оказаться более эффективным, если используются такие методы, как метод сопряжённых градиентов. Если имеется линейная зависимость столбцов матрицы J_r, метод итераций завершается неудачей, поскольку J_r^TJ_r становится вырожденной.

Пример

Вычисленная кривая с ${\hat {\beta }}_{1}=0.362$ и ${\hat {\beta }}_{2}=0.556$ (синим цветом) и наблюдаемые данные (красным цветом).

В этом примере используется алгоритм Гаусса — Ньютона для построения модели данных путём минимизации суммы квадратов отклонений данных и модели.

В экспериментальной биологии изучение связей между концентрацией субстрата [S] и скоростью реакции в реакции энзимомодуляции, были получены следующие данные.

i	1	2	3	4	5	6	7
[S]	0.038	0.194	0.425	0.626	1.253	2.500	3.740
скорость	0.050	0.127	0.094	0.2122	0.2729	0.2665	0.3317

Нужно найти кривую (функцию-модель) вида

скорость $={\frac {V_{\text{max}}[S]}{K_{M}+[S]}}$,

которая приближает наилучшим образом данные в смысле наименьших квадратов с параметрами $V_{\text{max}}$ и $K_{M}$, которые следует найти.

Обозначим через $x_{i}$ и $y_{i}$ значения [S] и скорость из таблицы, $i=1,\dots ,7$. Пусть $\beta _{1}=V_{\text{max}}$ и $\beta _{2}=K_{M}$. Мы будем искать $\beta _{1}$ и $\beta _{2}$, такие, что сумма квадратов отклонений

$r_{i}=y_{i}-{\frac {\beta _{1}x_{i}}{\beta _{2}+x_{i}}}\;(i=1,\dots ,7)$

минимальна.

Якобиан $\mathbf {J_{r}}$ вектора остатков $r_{i}$ по неизвестным $\beta _{j}$ — это $7\times 2$ матрица с $i$-ой строкой, имеющей элементы

${\frac {\partial r_{i}}{\partial \beta _{1}}}=-{\frac {x_{i}}{\beta _{2}+x_{i}}},\ {\frac {\partial r_{i}}{\partial \beta _{2}}}={\frac {\beta _{1}x_{i}}{\left(\beta _{2}+x_{i}\right)^{2}}}.$

Начав с начального приближения $\beta _{1}=0.9$ и $\beta _{2}=0.2$ после пяти итераций алгоритм Гаусса — Ньютона даёт оптимальные значения ${\hat {\beta }}_{1}=0.362$ и ${\hat {\beta }}_{2}=0.556$. Сумма квадратов остатков уменьшается от начального значения 1.445 до 0.00784 к пятой итерации. График справа показывает кривую с оптимальными параметрами.

Сходимость

Можно показать, что направление увеличения Δ является направлением спуска для S, и, если алгоритм сходится, пределом будет стационарная точка для S. Однако сходимость не гарантируется, даже когда начальная точка близка к решению, что происходит в методе Ньютона или методе BFGS при обычных условиях Фольфе.

Скорость сходимости алгоритма Гаусса — Ньютона близка к квадратичной. Алгоритм может сходиться медленнее или не сходиться совсем, если начальное приближение далеко от минимального или если матрица $\mathbf {J_{r}^{\mathsf {T}}J_{r}}$ плохо обусловлена. Например, представим задачу с $m=2$ уравнениями и $n=1$ переменной

${\begin{aligned}r_{1}(\beta )&=\beta +1\\r_{2}(\beta )&=\lambda \beta ^{2}+\beta -1.\end{aligned}}$

Полученное оптимальное решение — $\beta =0$. (Настоящий оптимум — $\beta =-1$ для $\lambda =2$, поскольку $S(0)=1^{2}+(-1)^{2}=2$, в то время как $S(-1)=0$.) Если $\lambda =0$, то задача, фактически, линейна и метод находит решение за одну итерацию. Если |λ| < 1, то метод сходится линейно и ошибка убывает со скоростью |λ| на каждой итерации. Однако, если |λ| > 1, то метод не сходится даже локально.

Алгоритм на основе метода Ньютона

Далее предполагается, что алгоритм Гаусса — Ньютона основан на методе Ньютона для минимизации функции с помощью аппроксимации. Как следствие, скорость сходимости алгоритма Гаусса — Ньютона может быть квадратична, если выполнены некоторые условия. В общем случае (при более слабых условиях), скорость сходимости может оказаться линейной.

Рекуррентное отношение метода Ньютона для минимизации функции S от параметров ${\boldsymbol {\beta }}$

${\boldsymbol {\beta }}^{(s+1)}={\boldsymbol {\beta }}^{(s)}-\mathbf {H} ^{-1}\mathbf {g} \,$

где g означает вектор-градиент функции S, а H обозначает гессиан функции S. Поскольку $S=\sum _{i=1}^{m}r_{i}^{2}$, градиент задаётся равенством

$g_{j}=2\sum _{i=1}^{m}r_{i}{\frac {\partial r_{i}}{\partial \beta _{j}}}.$

Элементы гессиана вычисляются путём дифференцирования элементов градиента $g_{j}$ по $\beta _{k}$

$H_{jk}=2\sum _{i=1}^{m}\left({\frac {\partial r_{i}}{\partial \beta _{j}}}{\frac {\partial r_{i}}{\partial \beta _{k}}}+r_{i}{\frac {\partial ^{2}r_{i}}{\partial \beta _{j}\partial \beta _{k}}}\right).$

Метод Гаусса — Ньютона получается путём отбрасывания второй производной (второго члена в выражении). То есть гессиан аппроксимируется

$H_{jk}\approx 2\sum _{i=1}^{m}J_{ij}J_{ik}$,

где $J_{ij}={\frac {\partial r_{i}}{\partial \beta _{j}}}$ — элементы якобиана J_r. Градиент и приближённый гессиан можно записать в матричной нотации

$\mathbf {g} =2\mathbf {J} _{\mathbf {r} }^{\mathsf {T}}\mathbf {r} ,\quad \mathbf {H} \approx 2\mathbf {J} _{\mathbf {r} }^{\mathsf {T}}\mathbf {J_{r}} .\,$

Эти выражения подставляются в рекуррентное отношение выше для получения операционных уравнений

${\boldsymbol {\beta }}^{(s+1)}={\boldsymbol {\beta }}^{(s)}+\Delta ;\quad \Delta =-\left(\mathbf {J_{r}} ^{\mathsf {T}}\mathbf {J_{r}} \right)^{-1}\mathbf {J_{r}} ^{\mathsf {T}}\mathbf {r} .$

Сходимость метода Гаусса — Ньютона в общем случае не гарантирована. Аппроксимация

$\left|r_{i}{\frac {\partial ^{2}r_{i}}{\partial \beta _{j}\partial \beta _{k}}}\right|\ll \left|{\frac {\partial r_{i}}{\partial \beta _{j}}}{\frac {\partial r_{i}}{\partial \beta _{k}}}\right|$,

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

1. Значения функции $r_{i}$ малы по величине, по меньшей мере, рядом с минимумом.
2. Функции лишь "слегка" нелинейны, то есть ${\frac {\partial ^{2}r_{i}}{\partial \beta _{j}\partial \beta _{k}}}$ относительно малы по величине.

Улучшенные версии

В методах Гаусса — Ньютона сумма квадратов остатков S может не уменьшаться на каждой итерации. Однако, поскольку Δ направлен в сторону уменьшения функции, если $S({\boldsymbol {\beta }}^{s})$ не является стационарной точкой, выполняется неравенство $S({\boldsymbol {\beta }}^{s}+\alpha \Delta )<S({\boldsymbol {\beta }}^{s})$ для достаточно малых $\alpha >0$. Таким образом, если обнаруживается расходимость, можно использовать долю $\alpha$ вектора приращения Δ в формуле обновления:

${\boldsymbol {\beta }}^{s+1}={\boldsymbol {\beta }}^{s}+\alpha \ \Delta$.

Другими словами, вектор приращения слишком длинный, но он указывает направление «спуска», так что если пройти лишь часть пути, можно уменьшить значение функции S. Оптимальное значение $\alpha$ может быть найдено с помощью алгоритма одномерного поиска, то есть величина $\alpha$ определяется путём нахождения значения, минимизирующего S с помощью одномерного поиска на интервале $0<\alpha <1$.

В случаях, когда в направлении вектора приращения оптимальная доля $\alpha$ близка к нулю, альтернативным методом отработки расходимости является использование алгоритма Левенберга — Марквардта, известного также как «метод доверительных областей». Нормальные уравнения, модифицированные таким образом, что вектор спуска поворачивается в направлении наискорейшего спуска,

$\left(\mathbf {J^{T}J+\lambda D} \right)\Delta =-\mathbf {J} ^{T}\mathbf {r}$,

где D — положительная диагональная матрица. Заметим, что если D является единичной матрицей E и $\lambda \to +\infty$, то $\lambda \Delta =\lambda \left(\mathbf {J^{E}J} +\lambda \mathbf {E} \right)^{-1}\left(-\mathbf {J} ^{T}\mathbf {r} \right)=\left(\mathbf {E} -\mathbf {J^{T}J} /\lambda +\cdots \right)\left(-\mathbf {J} ^{T}\mathbf {r} \right)\to -\mathbf {J} ^{T}\mathbf {r}$. Таким образом, направление Δ приближает направление отрицательного градиента $-\mathbf {J} ^{T}\mathbf {r}$.

Так называемый параметр Марквардта $\lambda$ может также быть оптимизирован путём линейного поиска, но смысла особого нет, поскольку вектор сдвига нужно каждый раз пересчитывать, когда меняется $\lambda$. Более эффективная стратегия такая. Если обнаруживается расхождение, увеличиваем параметр Марквардта, пока S убывает. Затем сохраняем значение между итерациями, но уменьшаем его, если возможно, пока не достигнем значения, когда параметр Марквардта не может быть обнулён. Минимизация S тогда становится стандартной минимизацией Гаусса — Ньютона.

Оптимизация больших задач

Для оптимизации большого размера метод Гаусса — Ньютона особенно интересен, поскольку часто (хотя, определённо, не всегда) матрица $\mathbf {J} _{\mathbf {r} }$ более разрежена, чем приближённый гессиан $\mathbf {J} _{\mathbf {r} }^{\mathsf {T}}\mathbf {J_{r}}$. В таких случаях шаг вычисления сам по себе, обычно, требует применения итеративного приближённого метода, такого как метод сопряжённых градиентов.

Чтобы этот подход работал, необходим как минимум эффективный метод вычисления произведения

$\mathbf {J} _{\mathbf {r} }^{\mathsf {T}}\mathbf {J_{r}} \mathbf {p}$

для некоторого вектора p. Для хранения разреженной матрицы практично хранить строки матрицы $\mathbf {J} _{\mathbf {r} }$ в сжатом виде (т.е. без нулевых элементов), что делает прямое вычисление приведённого выше произведения (ввиду транспозиции) сложным. Однако, если определить c_i как строку i матрицы $\mathbf {J} _{\mathbf {r} }$, выполняется следующее отношение:

$\mathbf {J} _{\mathbf {r} }^{\mathsf {T}}\mathbf {J_{r}} \mathbf {p} =\sum _{i}\mathbf {c} _{i}(\mathbf {c} _{i}\cdot \mathbf {p} )$

так что любая строка делает вклад в произведение аддитивно и независимо. Кроме того, это выражение хорошо изучено для применения параллельных вычислений. Заметим, что любая строка c_i является градиентом соответствующей невязки r_i. При учёте этого обстоятельства вышеприведённая формула подчёркивает факт, что невязки вносят вклад в результат независимо друг от друга.

Связанные алгоритмы

В квазиньютоновских методах, таких как методы Дэвидона, Флетчера и Пауэлла или Бройдена — Флетчера — Голдфарба — Шэнно (метод БФГШ) приближение полного гессиан ${\frac {\partial ^{2}S}{\partial \beta _{j}\partial \beta _{k}}}$ строится с помощью первых производных ${\frac {\partial r_{i}}{\partial \beta _{j}}}$ так, что после n уточнений метод по производительности близок к методу Ньютона. Заметим, что квазиньютоновские методы могут минимизировать вещественные функции общего вида, в то время как методы Гаусса — Ньютона, Левенберга — Марквардта и т.д. применимы только к нелинейным задачам наименьших квадратов.

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

Ссылки

Реализации

• Artelys Knitro. Система для решения нелинейных задач с импленментацией метода Гаусса — Ньютона. Система написана на C и имеет интерфейсы для C++/C#/Java/Python/MATLAB/R.