Русская Википедия:Алгоритм Гаусса — Ньютона

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

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

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

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

Описание

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

<math> S(\boldsymbol \beta)= \sum_{i=1}^m r_i^2(\boldsymbol \beta).</math>

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

<math> \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)}) </math>

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

<math> (\mathbf{J_r})_{ij} = \frac{\partial r_i (\boldsymbol \beta^{(s)})}{\partial \beta_j}</math>

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

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

<math> \boldsymbol \beta^{(s+1)} = \boldsymbol \beta^{(s)} - \left( \mathbf{J_r} \right)^{-1} \mathbf{r}(\boldsymbol \beta^{(s)}) </math>,

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

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

<math>r_i(\boldsymbol \beta)= y_i - f(x_i, \boldsymbol \beta).</math>

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

<math> \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)}). </math>

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

Замечания

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

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

<math>\mathbf{r}(\boldsymbol \beta)\approx \mathbf{r}(\boldsymbol \beta^s)+\mathbf{J_r}(\boldsymbol \beta^s)\Delta</math>,

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

<math>\mathbf{min}\|\mathbf{r}(\boldsymbol \beta^s)+\mathbf{J_r}(\boldsymbol \beta^s)\Delta\|_2^2</math>,

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

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

Пример

Файл:Gauss Newton illustration.png

Вычисленная кривая с <math>\hat\beta_1=0.362</math> и <math>\hat\beta_2=0.556</math> (синим цветом) и наблюдаемые данные (красным цветом).

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

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

Шаблон:Mvar	1	2	3	4	5	6	7
Шаблон:Math	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

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

скорость <math>=\frac{V_\text{max}[S]}{K_M+[S]}</math>,

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

Обозначим через <math>x_i</math> и <math>y_i</math> значения Шаблон:Math и скорость из таблицы, <math>i=1, \dots, 7</math>. Пусть <math>\beta_1=V_\text{max}</math> и <math>\beta_2=K_M</math>. Мы будем искать <math>\beta_1</math> и <math>\beta_2</math>, такие, что сумма квадратов отклонений

минимальна.

Якобиан <math>\mathbf{J_r}</math> вектора остатков <math>r_i</math> по неизвестным <math>\beta_j</math> — это <math>7\times 2</math> матрица с <math>i</math>-ой строкой, имеющей элементы

<math>\frac{\partial r_i}{\partial \beta_1}= -\frac{x_i}{\beta_2+x_i},\ \frac{\partial r_i}{\partial \beta_2}= \frac{\beta_1x_i}{\left(\beta_2+x_i\right)^2}.</math>

Начав с начального приближения <math>\beta_1=0.9</math> и <math>\beta_2=0.2</math> после пяти итераций алгоритм Гаусса — Ньютона даёт оптимальные значения <math>\hat\beta_1=0.362</math> и <math>\hat\beta_2=0.556</math>. Сумма квадратов остатков уменьшается от начального значения 1.445 до 0.00784 к пятой итерации. График справа показывает кривую с оптимальными параметрами.

Сходимость

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

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

<math> \begin{align}

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

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

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

Рекуррентное отношение метода Ньютона для минимизации функции S от параметров <math>\boldsymbol\beta</math>

<math> \boldsymbol\beta^{(s+1)} = \boldsymbol\beta^{(s)} - \mathbf H^{-1} \mathbf g \, </math>

где g означает вектор-градиент функции S, а H обозначает гессиан функции S. Поскольку <math> S = \sum_{i=1}^m r_i^2</math>, градиент задаётся равенством

<math>g_j=2\sum_{i=1}^m r_i\frac{\partial r_i}{\partial \beta_j}.</math>

Элементы гессиана вычисляются путём дифференцирования элементов градиента <math>g_j</math> по <math>\beta_k</math>

<math>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).</math>

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

<math>H_{jk}\approx 2\sum_{i=1}^m J_{ij}J_{ik}</math>,

где <math>J_{ij}=\frac{\partial r_i}{\partial \beta_j}</math> — элементы якобиана J_r. Градиент и приближённый гессиан можно записать в матричной нотации

<math>\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}.\,</math>

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

<math> \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}. </math>

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

<math>\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|</math>,

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

Значения функции <math>r_i</math> малы по величине, по меньшей мере, рядом с минимумом.
Функции лишь "слегка" нелинейны, то есть <math>\frac{\partial^2 r_i}{\partial \beta_j \partial \beta_k}</math> относительно малы по величине.

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

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

<math> \boldsymbol \beta^{s+1} = \boldsymbol \beta^s+\alpha\ \Delta</math>.

Другими словами, вектор приращения слишком длинный, но он указывает направление «спуска», так что если пройти лишь часть пути, можно уменьшить значение функции S. Оптимальное значение <math>\alpha</math> может быть найдено с помощью алгоритма [[|en]] (Line search), то есть величина <math>\alpha</math> определяется путём нахождения значения, минимизирующего S с помощью [[|en]] (Line search) на интервале <math>0<\alpha<1</math>.

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

<math>\left(\mathbf{J^TJ+\lambda D}\right)\Delta=-\mathbf{J}^T \mathbf{r}</math>,

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

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

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

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

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

<math>\mathbf{J}_\mathbf{r}^\mathsf{T}\mathbf{J_r} \mathbf p</math>

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

<math>\mathbf{J}_\mathbf{r}^\mathsf{T}\mathbf{J_r} \mathbf p = \sum_i \mathbf c_i (\mathbf c_i \cdot \mathbf p)</math>

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

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

В квазиньютоновских методах, таких как методы Шаблон:Не переведено 5 или Бройдена — Флетчера — Голдфарба — Шэнно (метод БФГШ) приближение полного гессиана <math>\frac{\partial^2 S}{\partial \beta_j \partial\beta_k}</math> строится с помощью первых производных <math>\frac{\partial r_i}{\partial\beta_j}</math> так, что после n уточнений метод по производительности близок к методу Ньютона. Заметим, что квазиньютоновские методы могут минимизировать вещественные функции общего вида, в то время как методы Гаусса — Ньютона, Левенберга — Марквардта и т.д. применимы только к нелинейным задачам наименьших квадратов.

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

Примечания

Шаблон:Примечания

Литература

Шаблон:Refbegin

Шаблон:Refend

Ссылки

Реализации

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

Шаблон:Методы оптимизации

Шаблон:Rq

Партнерские ресурсы
Криптовалюты	Обмен криптовалют - www.bestchange.ru Криптовалютная биржа CoinEx Криптовалютная биржа Binance HIVE OS - операционная система для майнинга e4pool - Мультивалютный пул для майнинга.
Магазины	AliExpress — глобальная виртуальная (в Интернете) торговая площадка, предоставляющая возможность покупать товары производителей из КНР; computeruniverse.net - Интернет-магазин компьютеров(Промо код 5 Евро на первую покупку:FWWC3ZKQ);
Хостинг	DigitalOcean - американский провайдер облачных инфраструктур, с главным офисом в Нью-Йорке и с центрами обработки данных по всему миру;
Разное	Викиум - Онлайн-тренажер для мозга Like Центр - Центр поддержки и развития предпринимательства. Gamersbay - лучший магазин по бустингу для World of Warcraft. Ноотропы OmniMind N°1 - Усиливает мозговую активность. Повышает мотивацию. Улучшает память. Санкт-Петербургская школа телевидения - это федеральная сеть образовательных центров, которая имеет филиалы в 37 городах России. Lingualeo.com — интерактивный онлайн-сервис для изучения и практики английского языка в увлекательной игровой форме. Junyschool (Джунискул) – международная школа программирования и дизайна для детей и подростков от 5 до 17 лет, где ученики осваивают компьютерную грамотность, развивают алгоритмическое и креативное мышление, изучают основы программирования и компьютерной графики, создают собственные проекты: игры, сайты, программы, приложения, анимации, 3D-модели, монтируют видео. Умназия - Интерактивные онлайн-курсы и тренажеры для развития мышления детей 6-13 лет SkillBox - это один из лидеров российского рынка онлайн-образования. Среди партнеров Skillbox ведущий разработчик сервисного дизайна AIC, медиа-компания Yoola, первое и самое крупное русскоязычное аналитическое агентство Tagline, онлайн-школа дизайна и иллюстрации Bang! Bang! Education, оператор PR-рынка PACO, студия рисования Draw&Go, агентство performance-маркетинга Ingate, scrum-студия Sibirix, имидж-лаборатория Персона. «Нетология» — это университет по подготовке и дополнительному обучению специалистов в области интернет-маркетинга, управления проектами и продуктами, дизайна, Data Science и разработки. В рамках Нетологии студенты получают ценные теоретические знания от лучших экспертов Рунета, выполняют практические задания на отработку полученных навыков, общаются с экспертами и единомышленниками. Познакомиться со всеми продуктами подробнее можно на сайте https://netology.ru, линейка курсов и профессий постоянно обновляется. StudyBay Brazil – это онлайн биржа для португалоговорящих студентов и авторов! Студент получает уникальную работу любого уровня сложности и больше свободного времени, в то время как у автора появляется дополнительный заработок и бесценный опыт. Автор24 — самая большая в России площадка по написанию учебных работ: контрольные и курсовые работы, дипломы, рефераты, решение задач, отчеты по практике, а так же любой другой вид работы. Сервис сотрудничает с более 70 000 авторов. Более 1 000 000 работ уже выполнено. StudyBay – это онлайн биржа для англоязычных студентов и авторов! Студент получает уникальную работу любого уровня сложности и больше свободного времени, в то время как у автора появляется дополнительный заработок и бесценный опыт.