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

Алгоритм Гаусса –Ньютона используется для решения нелинейных задач наименьших квадратов, что эквивалентно минимизации суммы квадратов значений функции. Это расширение метода Ньютона для нахождения минимума нелинейной функции . Поскольку сумма квадратов должна быть неотрицательной, алгоритм можно рассматривать как использующий метод Ньютона для итеративного приближения нулей компонентов суммы и, таким образом, минимизации суммы. В этом смысле алгоритм также является эффективным методом решения переопределенных систем уравнений. Его преимущество в том, что не требуются вторые производные, вычисление которых может быть сложным. ^[1]

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

Метод назван в честь математиков Карла Фридриха Гаусса и Исаака Ньютона и впервые появился в работе Гаусса 1809 года Theoria motus corporum coelestium insectionibus conicis solemambitum . ^[2]

Заданные функции (часто называемые остатками) переменных с алгоритмом Гаусса-Ньютона итеративно находят значение , которое минимизирует сумму квадратов ^[3]${\displaystyle м}$${\displaystyle {\textbf {r}}=(r_{1},\ldots ,r_{m})}$${\displaystyle n}$${\displaystyle {\boldsymbol {\beta }}=(\beta _{1},\ldots \beta _{n}),}$${\displaystyle m\geq n,}$${\displaystyle \бета}$$${\displaystyle S({\boldsymbol {\beta }})=\sum _{i=1}^{m}r_{i}({\boldsymbol {\beta }})^{2}.}$$

Начиная с первоначального предположения о минимуме, метод продолжается итерациями${\displaystyle {\boldsymbol {\beta }}^{(0)}}$$${\displaystyle {\boldsymbol {\beta }}^{(s+1)}={\boldsymbol {\beta }}^{(s)}-\left(\mathbf {J_{r}} ^{\operatorname {T} }\mathbf {J_{r}} \right)^{-1}\mathbf {J_{r}} ^{\operatorname {T} }\mathbf {r} \left({\boldsymbol {\beta }}^{(s)}\right),}$$

где, если r и β являются векторами-столбцами , элементы матрицы Якоби равны$${\displaystyle \left(\mathbf {J_{r}} \right)_{ij}={\frac {\partial r_{i}\left({\boldsymbol {\beta }}^{(s)}\right)}{\partial \beta _{j}}},}$$

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

На каждой итерации обновление можно найти, переставив предыдущее уравнение в следующие два шага:${\displaystyle \Delta ={\boldsymbol {\beta }}^{(s+1)}-{\boldsymbol {\beta }}^{(s)}}$

• ${\displaystyle \Delta =-\left(\mathbf {J_{r}} ^{\operatorname {T} }\mathbf {J_{r}} \right)^{-1}\mathbf {J_{r}} ^{\operatorname {T} }\mathbf {r} \left({\boldsymbol {\beta }}^{(s)}\right)}$
• ${\displaystyle \mathbf {J_{r}} ^{\operatorname {T} }\mathbf {J_{r}} \Delta =-\mathbf {J_{r}} ^{\operatorname {T} }\mathbf {r} \left({\boldsymbol {\beta }}^{(s)}\right)}$

При подстановках , , и это превращается в обычное матричное уравнение вида , которое затем можно решить различными методами (см. Примечания).${\textstyle A=\mathbf {J_{r}} ^{\operatorname {T} }\mathbf {J_{r}} }$${\displaystyle \mathbf {b} =-\mathbf {J_{r}} ^{\operatorname {T} }\mathbf {r} \left({\boldsymbol {\beta }}^{(s)}\right)}$${\displaystyle \mathbf {x} =\Delta }$${\displaystyle A\mathbf {x} =\mathbf {b} }$

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

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

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

При подгонке данных, когда целью является нахождение параметров таким образом, чтобы заданная модельная функция наилучшим образом соответствовала некоторым точкам данных , функции являются остатками :${\displaystyle {\boldsymbol {\beta }}}$${\displaystyle \mathbf {f} (\mathbf {x} ,{\boldsymbol {\beta }})}$${\displaystyle (x_{i},y_{i})}$${\displaystyle r_{i}}$$${\displaystyle r_{i}({\boldsymbol {\beta }})=y_{i}-f\left(x_{i},{\boldsymbol {\beta }}\right).}$$

Тогда метод Гаусса–Ньютона можно выразить через якобиан функции как${\displaystyle \mathbf {J_{f}} =-\mathbf {J_{r}} }$${\displaystyle \mathbf {f} }$$${\displaystyle {\boldsymbol {\beta }}^{(s+1)}={\boldsymbol {\beta }}^{(s)}+\left(\mathbf {J_{f}} ^{\operatorname {T} }\mathbf {J_{f}} \right)^{-1}\mathbf {J_{f}} ^{\operatorname {T} }\mathbf {r} \left({\boldsymbol {\beta }}^{(s)}\right).}$$

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

Предположение m ≥ n в формулировке алгоритма необходимо, так как в противном случае матрица необратима и нормальные уравнения не могут быть решены (по крайней мере, однозначно).${\displaystyle \mathbf {J_{r}} ^{T}\mathbf {J_{r}} }$

Алгоритм Гаусса–Ньютона может быть получен путем линейной аппроксимации вектора функций r _i . Используя теорему Тейлора , мы можем записать на каждой итерации:$${\displaystyle \mathbf {r} ({\boldsymbol {\beta }})\approx \mathbf {r} \left({\boldsymbol {\beta }}^{(s)}\right)+\mathbf {J_{r}} \left({\boldsymbol {\beta }}^{(s)}\right)\Delta }$$

с . Задача нахождения минимизации суммы квадратов правой части; т.е.${\displaystyle \Delta ={\boldsymbol {\beta }}-{\boldsymbol {\beta }}^{(s)}}$${\displaystyle \Delta }$$${\displaystyle \min \left\|\mathbf {r} \left({\boldsymbol {\beta }}^{(s)}\right)+\mathbf {J_{r}} \left({\boldsymbol {\beta }}^{(s)}\right)\Delta \right\|_{2}^{2},}$$

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

Нормальные уравнения — это n одновременных линейных уравнений относительно неизвестных приращений . Их можно решить за один шаг, используя разложение Холецкого или, лучше, QR-факторизацию . Для больших систем итерационный метод , такой как метод сопряженных градиентов , может быть более эффективным. Если между столбцами J _r есть линейная зависимость , итерации не будут выполнены, так как становится сингулярным.${\displaystyle \Delta }$${\displaystyle \mathbf {J_{r}} }$${\displaystyle \mathbf {J_{r}} ^{T}\mathbf {J_{r}} }$

Если является сложным, следует использовать сопряженную форму: .${\displaystyle \mathbf {r} }$${\displaystyle \mathbf {r} :\mathbb {C} ^{n}\to \mathbb {C} }$${\displaystyle \left({\overline {\mathbf {J_{r}} }}^{\operatorname {T} }\mathbf {J_{r}} \right)^{-1}{\overline {\mathbf {J_{r}} }}^{\operatorname {T} }}$

В этом примере алгоритм Гаусса–Ньютона будет использоваться для подгонки модели к некоторым данным путем минимизации суммы квадратов ошибок между данными и прогнозами модели.

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

Требуется найти кривую (модельную функцию) вида$${\displaystyle {\text{rate}}={\frac {V_{\text{max}}\cdot [S]}{K_{M}+[S]}}}$$

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

Обозначим через и значения [ S ] и скорости соответственно, причем . Пусть и . Найдем и такие, что сумма квадратов остатков${\displaystyle x_{i}}$${\displaystyle y_{i}}$${\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}}},\quad (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}}};\quad {\frac {\partial r_{i}}{\partial \beta _{2}}}={\frac {\beta _{1}\cdot x_{i}}{\left(\beta _{2}+x_{i}\right)^{2}}}.}$$

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

Итерация Гаусса-Ньютона гарантированно сходится к точке локального минимума при выполнении 4 условий: ^[4] Функции дважды непрерывно дифференцируемы в открытом выпуклом множестве , якобиан имеет полный ранг столбца, начальная итерация близка к , а значение локального минимума мало. Сходимость квадратичная, если .${\displaystyle {\hat {\beta }}}$${\displaystyle r_{1},\ldots ,r_{m}}$${\displaystyle D\ni {\hat {\beta }}}$${\displaystyle \mathbf {J} _{\mathbf {r} }({\hat {\beta }})}$${\displaystyle \beta ^{(0)}}$${\displaystyle {\hat {\beta }}}$${\displaystyle |S({\hat {\beta }})|}$${\displaystyle |S({\hat {\beta }})|=0}$

Можно показать ^[5] , что приращение Δ является направлением спуска для S , и, если алгоритм сходится, то пределом является стационарная точка S. Однако для большого минимального значения сходимость не гарантируется, даже локальная сходимость , как в методе Ньютона , или сходимость при обычных условиях Вульфа. ^[6]${\displaystyle |S({\hat {\beta }})|}$

Скорость сходимости алгоритма Гаусса–Ньютона может приближаться к квадратичной . ^[7] Алгоритм может сходиться медленно или не сходиться вообще, если начальное предположение далеко от минимума или матрица плохо обусловлена . Например, рассмотрим задачу с уравнениями и переменной, заданной как${\displaystyle \mathbf {J_{r}^{\operatorname {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}}}$$

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

Итерация Гаусса-Ньютона является эффективным методом решения переопределенных систем уравнений в виде с и где — обратная матрица Мура-Пенроуза (также известная как псевдообратная матрица ) матрицы Якоби . Ее можно считать расширением метода Ньютона , и она обладает той же локальной квадратичной сходимостью ^[4] к изолированным регулярным решениям. $${\displaystyle \mathbf {x} ^{(k+1)}=\mathbf {x} ^{(k)}-J(\mathbf {x} ^{(k)})^{\dagger }\mathbf {f} (\mathbf {x} ^{(k)})\,,\quad k=0,1,\ldots }$$${\displaystyle \mathbf {f} (\mathbf {x} )=\mathbf {0} }$$${\displaystyle \mathbf {f} (\mathbf {x} )={\begin{bmatrix}f_{1}(x_{1},\ldots ,x_{n})\\\vdots \\f_{m}(x_{1},\ldots ,x_{n})\end{bmatrix}}}$$${\displaystyle m>n}$${\displaystyle J(\mathbf {x} )^{\dagger }}$ ${\displaystyle J(\mathbf {x} )}$${\displaystyle \mathbf {f} (\mathbf {x} )}$

Если решение не существует, но начальная итерация близка к точке , в которой сумма квадратов достигает небольшого локального минимума, итерация Гаусса-Ньютона линейно сходится к . Точку часто называют решением наименьших квадратов переопределенной системы.${\displaystyle \mathbf {x} ^{(0)}}$${\displaystyle {\hat {\mathbf {x} }}=({\hat {x}}_{1},\ldots ,{\hat {x}}_{n})}$${\textstyle \sum _{i=1}^{m}|f_{i}(x_{1},\ldots ,x_{n})|^{2}\equiv \|\mathbf {f} (\mathbf {x} )\|_{2}^{2}}$${\displaystyle {\hat {\mathbf {x} }}}$${\displaystyle {\hat {\mathbf {x} }}}$

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

Рекуррентное соотношение для метода Ньютона для минимизации функции S параметров имеет вид${\displaystyle {\boldsymbol {\beta }}}$$${\displaystyle {\boldsymbol {\beta }}^{(s+1)}={\boldsymbol {\beta }}^{(s)}-\mathbf {H} ^{-1}\mathbf {g} ,}$$

где g обозначает вектор градиента S , а H обозначает матрицу Гессе S.

Так как , градиент задается выражением${\textstyle 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},}$$

где — элементы якобиана J _r . Обратите внимание, что когда точный гессиан оценивается вблизи точного соответствия, мы имеем почти ноль , поэтому второй член также становится близким к нулю, что оправдывает приближение. Градиент и приближенный гессиан можно записать в матричной записи как${\textstyle J_{ij}={\partial r_{i}}/{\partial \beta _{j}}}$${\displaystyle r_{i}}$$${\displaystyle \mathbf {g} =2{\mathbf {J} _{\mathbf {r} }}^{\operatorname {T} }\mathbf {r} ,\quad \mathbf {H} \approx 2{\mathbf {J} _{\mathbf {r} }}^{\operatorname {T} }\mathbf {J_{r}} .}$$

Эти выражения подставляются в рекуррентное соотношение выше для получения рабочих уравнений$${\displaystyle {\boldsymbol {\beta }}^{(s+1)}={\boldsymbol {\beta }}^{(s)}+\Delta ;\quad \Delta =-\left(\mathbf {J_{r}} ^{\operatorname {T} }\mathbf {J_{r}} \right)^{-1}\mathbf {J_{r}} ^{\operatorname {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|}$$

которое необходимо для того, чтобы можно было игнорировать производные второго порядка, может быть справедливо в двух случаях, для которых следует ожидать сходимости: ^[10]

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

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

Другими словами, вектор приращения слишком длинный, но он все еще указывает «вниз», поэтому прохождение только части пути уменьшит целевую функцию S . Оптимальное значение для можно найти с помощью алгоритма линейного поиска , то есть величина определяется путем нахождения значения, которое минимизирует S , обычно с использованием метода прямого поиска в интервале или линейного поиска с возвратом, такого как поиск Armijo-line . Как правило, следует выбирать так, чтобы он удовлетворял условиям Вульфа или условиям Голдштейна . ^[11]${\displaystyle \alpha }$${\displaystyle \alpha }$${\displaystyle 0<\alpha <1}$${\displaystyle \alpha }$

В случаях, когда направление вектора сдвига таково, что оптимальная доля α близка к нулю, альтернативным методом обработки расхождения является использование алгоритма Левенберга–Марквардта , метода доверительной области . ^[3] Нормальные уравнения модифицируются таким образом, что вектор приращения поворачивается в направлении наискорейшего спуска ,$${\displaystyle \left(\mathbf {J^{\operatorname {T} }J+\lambda D} \right)\Delta =-\mathbf {J} ^{\operatorname {T} }\mathbf {r} ,}$$

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

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

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

Чтобы этот подход заработал, нужен как минимум эффективный метод вычисления произведения$${\displaystyle {\mathbf {J} _{\mathbf {r} }}^{\operatorname {T} }\mathbf {J_{r}} \mathbf {p} }$$

для некоторого вектора p . При хранении разреженных матриц обычно практично хранить строки в сжатом виде (например, без нулевых записей), что делает прямое вычисление вышеуказанного произведения сложным из-за транспонирования. Однако, если определить c _i как строку i матрицы , выполняется следующее простое соотношение:${\displaystyle \mathbf {J} _{\mathbf {r} }}$${\displaystyle \mathbf {J} _{\mathbf {r} }}$$${\displaystyle {\mathbf {J} _{\mathbf {r} }}^{\operatorname {T} }\mathbf {J_{r}} \mathbf {p} =\sum _{i}\mathbf {c} _{i}\left(\mathbf {c} _{i}\cdot \mathbf {p} \right),}$$

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

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

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

Следующая реализация в Julia предоставляет один метод, который использует предоставленный якобиан, и другой метод, вычисляющий с автоматическим дифференцированием .

"""  гауссньютон(r,J,β₀,макситер,тол)Выполнить оптимизацию Гаусса-Ньютона для минимизации остаточной функции `r` с якобианом `J`, начиная с `β₀`. Алгоритм завершается, когда норма шага меньше `tol` или после итераций `maxiter`. """ function gaussnewton ( r , J , β₀ , maxiter , tol ) β = copy ( β₀ ) for _ in 1 : maxiter Jβ = J ( β ); Δ = - ( Jβ '* Jβ ) \ ( Jβ '* r ( β )) β += Δ if sqrt ( sum ( abs2 , Δ )) < tol break end end return β end                             import AbstractDifferentiation as AD , Zygote backend = AD.ZygoteBackend () # доступны другие бэкенды       """  гауссньютон(r,β₀,макситер,тол)Выполнить оптимизацию Гаусса-Ньютона для минимизации остаточной функции `r`, начиная с `β₀`. Соответствующий якобиан вычисляется с помощью автоматического дифференцирования. Алгоритм завершается, когда норма шага становится меньше `tol` или после итераций `maxiter`. """ function gaussnewton ( r , β₀ , maxiter , tol ) β = copy ( β₀ ) for _ in 1 : maxiter rβ , Jβ = AD . value_and_jacobian ( backend , r , β ) Δ = - ( Jβ [ 1 ] '* Jβ [ 1 ]) \ ( Jβ [ 1 ] '* rβ ) β += Δ if sqrt ( sum ( abs2 , Δ )) < tol break end end return β end                              

• Бьёрк, А. (1996). Численные методы решения задач наименьших квадратов . SIAM, Филадельфия. ISBN 0-89871-360-9.
• Флетчер, Роджер (1987). Практические методы оптимизации (2-е изд.). Нью-Йорк: John Wiley & Sons . ISBN 978-0-471-91547-8..
• Нокедаль, Хорхе; Райт, Стивен (1999). Численная оптимизация . Нью-Йорк: Springer. ISBN 0-387-98793-2.

• Вероятность, статистика и оценка. Алгоритм подробно описан и применен к биологическому эксперименту, обсуждаемому в качестве примера в этой статье (стр. 84 с неопределенностями оценочных значений).

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