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

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

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

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

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

Описание

Данный м функции р = (р₁, …, р_м) (часто называемые остатками) п переменные β = (β₁, …, β_п), с м ≥ п, алгоритм Гаусса – Ньютона итеративно находит значение переменных, которое минимизирует сумму квадратов^[3]

${ 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} left ({ boldsymbol { beta }} ^ {(s)} right),}$

где, если р и β находятся вектор-столбец, записи Матрица якобиана находятся

${ displaystyle left ( mathbf {J_ {r}} right) _ {ij} = { frac { partial r_ {i} left ({ boldsymbol { beta}} ^ {(s)} справа)} { partial beta _ {j}}},}$

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

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

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

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

При подборе данных, цель которого - найти параметры β так что данная модельная функция у = ж(Икс, β) лучше всего подходит для некоторых точек данных (Икс_я, у_я) функции р_я являются остатки:

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

Тогда метод Гаусса – Ньютона можно выразить через якобиан J_ж функции ж в качестве

${ 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} left ({ boldsymbol { beta }} ^ {(s)} right).}$

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

Примечания

Предположение м ≥ п в формулировке алгоритма необходимо, иначе матрица J_р^ТJ_р не обратима, и нормальные уравнения не могут быть решены (по крайней мере, однозначно).

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

${ displaystyle mathbf {r} ({ boldsymbol { beta}}) приблизительно mathbf {r} left ({ boldsymbol { beta}} ^ {(s)} right) + mathbf {J_ {r}} left ({ boldsymbol { beta}} ^ {(s)} right) Delta}$

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

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

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

Нормальные уравнения: п совместных линейных уравнений с неизвестными приращениями Δ. Их можно решить за один шаг, используя Разложение Холецкого, или, лучше, QR-факторизация из J_р. Для больших систем итерационный метод, такой как сопряженный градиент метод, может быть более эффективным. Если существует линейная зависимость между столбцами J_р, итерации не удастся, так как J_р^ТJ_р становится единичным.

Когда ${ displaystyle mathbf {r}}$ сложно ${ displaystyle mathbf {r}}$:C^п${ displaystyle rightarrow}$C следует использовать конъюгированную форму: ${ displaystyle left ({ overline { mathbf {J_ {f}}}} ^ { mathsf {T}} mathbf {J_ {f}} right) ^ {- 1} { overline { mathbf {J_ {f}}}} ^ { mathsf {T}}}$.

Пример

Расчетная кривая, полученная с ${ displaystyle { hat { beta}} _ {1} = 0,362}$ и ${ displaystyle { hat { beta}} _ {2} = 0,556}$ (синим цветом) по сравнению с наблюдаемыми данными (красным).

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

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

я	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 { text {rate}} = { frac {V _ { text {max}} [S]} {K_ {M} + [S]}}}$

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

Обозначим через ${ displaystyle x_ {i}}$ и ${ displaystyle y_ {i}}$ значение [S] и ставка из таблицы, ${ Displaystyle я = 1, точки, 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, точки, 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 после пятой итерации. График на рисунке справа показывает кривую, определенную моделью для оптимальных параметров с наблюдаемыми данными.

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

Это можно показать^[4] что приращение Δ является направление спуска за S, и, если алгоритм сходится, то предел равен стационарная точка из S. Однако сходимость не гарантируется, даже локальная конвергенция как в Метод Ньютона, или сходимость в обычных условиях Вульфа.^[5]

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

${ displaystyle { begin {align} r_ {1} ( beta) & = beta +1, r_ {2} ( beta) & = lambda beta ^ {2} + beta -1. конец {выровнено}}}$

Оптимум на ${ 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, то метод даже не сходится локально.^[7]

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

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

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

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

куда грамм обозначает вектор градиента из S, и ЧАС обозначает Матрица Гессе из S.

С ${ Displaystyle S = сумма _ {я = 1} ^ {m} r_ {я} ^ {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} приблизительно 2 sum _ {i = 1} ^ {m} J_ {ij} J_ {ik},}$

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

${ displaystyle mathbf {g} = 2 { mathbf {J} _ { mathbf {r}}} ^ { mathsf {T}} mathbf {r}, quad mathbf {H} приблизительно 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 |}$

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

1. Значения функции ${ displaystyle r_ {i}}$ невелики по величине, по крайней мере, около минимума.
2. Функции только "умеренно" нелинейны, так что ${ displaystyle { 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)$ для всех достаточно малых ${ displaystyle alpha> 0}$. Таким образом, если происходит расхождение, одним из решений является использование дроби ${ displaystyle alpha}$ вектора приращения Δ в формуле обновления:

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

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

В случаях, когда направление вектора сдвига таково, что оптимальная доля α близка к нулю, альтернативным методом обработки дивергенции является использование Алгоритм Левенберга – Марквардта, а регион доверия метод.^[3] Нормальные уравнения изменяются таким образом, что вектор приращения поворачивается в направлении крутой спуск,

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

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

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

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