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

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

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

Описание

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

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

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

${\displaystyle {\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 и β как вектор-столбцы, элементы матрицы Якоби равны

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

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

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

${\displaystyle {\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 являются

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

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

${\displaystyle {\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)}).}$

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

Замечания

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

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

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

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

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

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

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

Пример

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

В экспериментальной биологии изучение связей между концентрацией субстрата [ 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

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

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

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

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

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

минимальна.

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

${\displaystyle {\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}}}.}$

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

Сходимость

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

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

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

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

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

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

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

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

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

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

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

${\displaystyle 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).}$

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

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

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

${\displaystyle \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}} .\,}$

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

${\displaystyle {\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} .}$

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

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

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

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

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

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

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

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

где D — положительная диагональная матрица. Заметим, что если D является единичной матрицей E и ${\displaystyle \lambda \to +\infty }$ , то ${\displaystyle \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} }$ . Таким образом, направление Δ приближает направление отрицательного градиента ${\displaystyle -\mathbf {J} ^{T}\mathbf {r} }$ .

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

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

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

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

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

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

${\displaystyle \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 . При учёте этого обстоятельства вышеприведённая формула подчёркивает факт, что невязки вносят вклад в результат независимо друг от друга.

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

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

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

Примечания

Литература

Ссылки

Реализации

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