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

Аппроксимация зашумленной кривой моделью асимметричного пика $f_{\beta }(x)$ с параметрами $\beta$ путем минимизации суммы квадратов остатков $r_{i}(\beta )=y_{i}-f_{\beta }(x_{i})$ в точках сетки $x_{i}$, используя алгоритм Гаусса – Ньютона.
Вверху: необработанные данные и модель.
Внизу: Эволюция нормализованной суммы квадратов ошибок.

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

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

Метод назван в честь математиков Карла Фридриха Гаусса и Исаака Ньютона и впервые появился в работе Гаусса 1809 года «Теория движения небесных тел в конических сечениях, окружающих Солнце» .

Описание

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

Начиная с первоначального предположения ${\boldsymbol {\beta }}^{(0)}$ для минимума метод выполняется итерациями $${\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 и β — векторы-столбцы , элементы матрицы Якоби равны $$\left(\mathbf {J_{r}} \right)_{ij}={\frac {\partial r_{i}\left({\boldsymbol {\beta }}^{(s)}\right)}{\partial \beta _{j}}},$$

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

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

• $\Delta =-\left(\mathbf {J_{r}} ^{\operatorname {T} }\mathbf {J_{r}} \right)^{-1}\mathbf {J_{r}} ^{\operatorname {T} }\mathbf {r} \left({\boldsymbol {\beta }}^{(s)}\right)$
• $\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}} }$, $\mathbf {b} =-\mathbf {J_{r}} ^{\operatorname {T} }\mathbf {r} \left({\boldsymbol {\beta }}^{(s)}\right)$, и $\mathbf {x} =\Delta$, это превращается в обычное матричное уравнение вида $A\mathbf {x} =\mathbf {b}$, которую затем можно решить различными методами (см. Примечания ).

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

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

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

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

Тогда метод Гаусса – Ньютона можно выразить через якобиан $\mathbf {J_{f}} =-\mathbf {J_{r}}$ функции $\mathbf {f}$ как $${\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).$$

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

Примечания

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

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

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

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

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

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

Пример

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

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

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

я	1	2	3	4	5	6	7
[ С ]	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

Требуется найти кривую (модельную функцию) вида $${\text{rate}}={\frac {V_{\text{max}}\cdot [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}}},\quad (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}}};\quad {\frac {\partial r_{i}}{\partial \beta _{2}}}={\frac {\beta _{1}\cdot 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 после пятой итерации. График на рисунке справа показывает кривую, определенную моделью для оптимальных параметров с учетом наблюдаемых данных.

Свойства сходимости

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

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

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

Решение переопределенных систем уравнений

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

Если решение не существует, но существует начальная итерация $\mathbf {x} ^{(0)}$ находится рядом с точкой ${\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}}$ достигает небольшого локального минимума, итерация Гаусса-Ньютона линейно сходится к ${\hat {\mathbf {x} }}$. Суть ${\hat {\mathbf {x} }}$ часто называют наименьших квадратов решением переопределенной системы методом .

Вывод из метода Ньютона

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

Рекуррентное соотношение для метода Ньютона минимизации функции S параметров ${\boldsymbol {\beta }}$ является $${\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}}$, градиент определяется выражением $$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},$$

где ${\textstyle J_{ij}={\partial r_{i}}/{\partial \beta _{j}}}$ являются элементами якобиана J _r . Обратите внимание, что когда точный гессиан оценивается вблизи точного соответствия, мы имеем почти нулевое значение. $r_{i}$, поэтому второй член также становится близким к нулю, что оправдывает аппроксимацию. Градиент и приближенный гессиан можно записать в матричной записи как $$\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}} .$$

Эти выражения подставляются в приведенное выше рекуррентное соотношение для получения рабочих уравнений $${\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} .$$

Сходимость метода Гаусса–Ньютона не гарантируется во всех случаях. Приближение $$\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. Функции являются лишь «слегка» нелинейными, так что ${\textstyle {\frac {\partial ^{2}r_{i}}{\partial \beta _{j}\partial \beta _{k}}}}$ имеет относительно небольшую величину.

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

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

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

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

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

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

Масштабная оптимизация

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

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

для некоторого вектора p . При хранении разреженной матрицы , как правило, практично хранить строки $\mathbf {J} _{\mathbf {r} }$ в сжатой форме (например, без нулевых записей), что затрудняет прямое вычисление вышеуказанного произведения из-за транспонирования. Однако если определить c _i как строку i матрицы $\mathbf {J} _{\mathbf {r} }$, имеет место следующее простое соотношение: $${\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 ), оценка полного гессиана ${\textstyle {\frac {\partial ^{2}S}{\partial \beta _{j}\partial \beta _{k}}}}$ строится численно с использованием первых производных ${\textstyle {\frac {\partial r_{i}}{\partial \beta _{j}}}}$ только для того, чтобы после n циклов уточнения метод по производительности максимально приблизился к методу Ньютона. Обратите внимание, что квазиньютоновские методы могут минимизировать общие вещественные функции, тогда как методы Гаусса – Ньютона, Левенберга – Марквардта и т. Д. подходят только для нелинейных задач наименьших квадратов.

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

Пример реализации

Юлия

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

"""
    gaussnewton(r,J,β₀,maxiter,tol)

Perform Gauss-Newton optimization to minimize the residual function `r` with Jacobian `J` starting from `β₀`. The algorithm terminates when the norm of the step is less than `tol` or after `maxiter` iterations.
"""
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() # other backends are available

"""
    gaussnewton(r,β₀,maxiter,tol)

Perform Gauss-Newton optimization to minimize the residual function `r` starting from `β₀`. The relevant Jacobian is calculated using automatic differentiation. The algorithm terminates when the norm of the step is less than `tol` or after `maxiter` iterations.
"""
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). Численные методы решения задач наименьших квадратов . СИАМ, Филадельфия. ISBN  0-89871-360-9 .
• Флетчер, Роджер (1987). Практические методы оптимизации (2-е изд.). Нью-Йорк: Джон Уайли и сыновья . ISBN  978-0-471-91547-8 . .
• Носедаль, Хорхе; Райт, Стивен (1999). Численная оптимизация . Нью-Йорк: Спрингер. ISBN  0-387-98793-2 .

Внешние ссылки

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

Реализации

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