Наименьшие квадраты (подбор кривой модели) алгоритмы

Определение наименьших квадратов

Наименьшие квадраты, в целом, являются задачей нахождения векторный x, который является локальным минимизатором к функции, которая является суммой квадратов, возможно подвергните некоторым ограничениям:

$\min_{x} {‖ F (x) ‖}_{2}^{2} = \min_{x} \sum_{i} F_{i}^{2} (x)$

таким образом, что A·x ≤ b, Aeq·x = beq, lb ≤ x ≤ ub.

Существует несколько решателей Optimization Toolbox™, доступных для различных типов F (x) и различные типы ограничений:

Решатель	F (x)	Ограничения
`mldivide`	C· – d	'none'
`lsqnonneg`	C· – d	x ≥ 0
`lsqlin`	C· – d	Связанный, линейный
`lsqnonlin`	Общий F (x)	Связанный
`lsqcurvefit`	F (x, xdata) – ydata	Связанный

Существуют алгоритмы на пять наименьших квадратов в решателях Optimization Toolbox, в дополнение к алгоритмам, используемым в mldivide:

lsqlin внутренняя точка
lsqlin активный набор
"Доверительная область, отражающая" (нелинейный или линейный метод наименьших квадратов)
Levenberg-Marquardt (нелинейный метод наименьших квадратов)
Алгоритм, используемый lsqnonneg

Все алгоритмы кроме lsqlin активный набор является крупномасштабным; смотрите Крупномасштабный по сравнению с Алгоритмами Средней шкалы. Для общего обзора методов нелинейного метода наименьших квадратов смотрите Денниса [8]. Определенные детали о методе Levenberg-Marquardt могут быть найдены в Moré [28].

Линейный метод наименьших квадратов: внутренняя точка или Активный Набор

lsqlin 'interior-point' алгоритм использует выпуклый внутренней точкой quadprog Алгоритм, и lsqlin 'active-set' алгоритм использует активный набор quadprog алгоритм. quadprog описание задачи должно минимизировать квадратичную функцию

$\min_{x} \frac{1}{2} x^{T} H x + c^{T} x$

подвергните линейным ограничениям и связанным ограничениям. lsqlin функция минимизирует квадратичную норму векторного Cx – d, удовлетворяющего линейным ограничениям и связанным ограничениям. Другими словами, lsqlin минимизирует

$\begin{matrix} {‖ C x - d ‖}_{2}^{2} = {(C x - d)}^{T} (C x - d) \\ = (x^{T} C^{T} - d^{T}) (C x - d) \\ = (x^{T} C^{T} C x) - (x^{T} C^{T} d - d^{T} C x) + d^{T} d \\ = \frac{1}{2} x^{T} (2 C^{T} C) x + {(- 2 C^{T} d)}^{T} x + d^{T} d . \end{matrix}$

Это помещается quadprog среда путем установки матрицы H на 2^CTC и вектора c к ^(–2CTd). (Вектор свободных членов ^dTd не оказывает влияния на местоположение минимума.) После этой переформулировки lsqlin проблема, quadprog вычисляет решение.

Примечание

quadprog 'interior-point-convex' алгоритм имеет два пути выполнения кода. Это берет тот, когда матрица Гессиана, из которой H является обычной (полной) матрицей, удваивается, и это берет другой, когда H является разреженной матрицей. Для получения дополнительной информации разреженного типа данных, смотрите Разреженные матрицы. Обычно алгоритм быстрее для больших проблем, которые имеют относительно немного ненулевых условий, когда вы задаете H как sparse. Точно так же алгоритм быстрее для небольших или относительно плотных проблем, когда вы задаете H как full.

Доверительная область отражающие наименьшие квадраты

Доверительная область отражающий алгоритм наименьших квадратов

Многие методы, используемые в решателях Optimization Toolbox, основаны на доверительных областях, простой все же мощной концепции в оптимизации.

Чтобы изучить подход доверительной области к оптимизации, рассмотрите задачу безусловной минимизации, минимизируйте f (x), где функция берет аргументы вектора и возвращает скаляры. Предположим, что вы - в точке x в n - пробел, и вы хотите улучшиться, т.е. переместиться в точку с более низким значением функции. Основная идея состоит в том, чтобы аппроксимировать f более простым функциональным q, который обоснованно отражает поведение функционального f в окружении N вокруг точки x. Это окружение является доверительной областью. Испытательный шаг s вычисляется путем минимизации (или приблизительно минимизации) по N. Это - подпроблема доверительной области,

\min_{s} {q (s), s \in N} .

(1)

Текущая точка обновляется, чтобы быть x + s если f (x + s) < f (x); в противном случае текущая точка остается неизменной, и N, область доверия, уменьшается, и испытательный расчет шага повторяется.

Ключевые вопросы в определении определенного подхода доверительной области к минимизации f (x) состоят в том, как выбрать и вычислить приближение q (заданный в текущей точке x), как выбрать и изменить доверительную область N, и как точно решить подпроблему доверительной области. Этот раздел фокусируется на неограниченной проблеме. Более поздние разделы обсуждают дополнительные осложнения из-за присутствия ограничений на переменные.

В стандартном методе доверительной области ([48]), квадратичное приближение q задан первыми двумя сроками приближения Тейлора к F в x; окружение N является обычно сферическим или эллипсоидальным в форме. Математически подпроблема доверительной области обычно утверждается

\min {\frac{1}{2} s^{T} H s + s^{T} g таким образом , что ‖ D s ‖ \leq Δ},

(2)

где g является градиентом f в текущей точке x, H является матрицей Гессиана (симметрическая матрица вторых производных), D является диагональной матрицей масштабирования, Δ является положительной скалярной величиной и ∥. ∥ является 2-нормой. Хорошие алгоритмы существуют для решения уравнения 2 (см. [48]); такие алгоритмы обычно включают расчет всех собственных значений H, и процесс Ньютона применился к характеристическому уравнению

$\frac{1}{Δ} - \frac{1}{‖ s ‖} = 0.$

Такие алгоритмы предоставляют точное решение уравнения 2. Однако они требуют времени, пропорционального нескольким факторизациям H. Поэтому для проблем доверительной области другой подход необходим. Несколько приближений и эвристических стратегий, на основе уравнения 2, были предложены в литературе ([42] и [50]). Подход приближения, сопровождаемый в решателях Optimization Toolbox, должен ограничить подпроблему доверительной области двумерным подпространством S ([39] и [42]). Однажды подпространство S был вычислен, работа, чтобы решить уравнение 2 тривиальна, даже если полная информация о собственном значении/собственном векторе необходима (поскольку в подпространстве, проблема только двумерна). Доминирующая работа теперь переключила к определению подпространства.

Двумерное подпространство S определяется при помощи предобусловленного процесса метода сопряженных градиентов, описанного ниже. Решатель задает S как линейный пробел, заполненный s ₁ и s ₂, где s ₁ является в направлении градиента g, и s ₂ является любой аппроксимированным направлением Ньютона, i.e., решение

H \cdot s_{2} = - g,

(3)

или направление отрицательной кривизны,

s_{2}^{T} \cdot H \cdot s_{2} < 0.

(4)

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

Эскиз безусловной минимизации с помощью идей доверительной области теперь легко дать:

Сформулируйте двумерную подпроблему доверительной области.
Решите уравнение 2, чтобы определить испытательный шаг s.
Если f (x + s) <f (x), то x = x + s.
Настройте Δ.

Эти четыре шага повторяются до сходимости. Размерность доверительной области Δ настроена согласно стандартным правилам. В частности, это уменьшено, если испытательный шаг не принят, т.е. f (x + s) ≥ f (x). См. [46] и [49] для обсуждения этого аспекта.

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

Крупномасштабный нелинейный метод наименьших квадратов

Важный особый случай для f (x) является проблемой нелинейного метода наименьших квадратов

\min_{x} \sum_{i} f_{i}^{2} (x) = \min_{x} {‖ F (x) ‖}_{2}^{2},

(5)

где F (x) является функцией с векторным знаком с i компонента F (x), равный _fi (x). Основной метод, используемый, чтобы решить эту задачу, эквивалентен в общем случае, описанном в Методах Доверительной области для Нелинейной Минимизации. Однако структура проблемы нелинейного метода наименьших квадратов используется, чтобы улучшить КПД. В частности, аппроксимированное направление Ньютона Гаусса, т.е. решение s к

\min {‖ J s + F ‖}_{2}^{2},

(6)

(где J является якобианом F (x)), используется, чтобы помочь задать двумерное подпространство S. Вторые производные функционального _fi компонента (x) не используются.

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

$J^{T} J s = - J^{T} F,$

несмотря на то, что нормальные уравнения явным образом не формируются.

Крупномасштабный линейный метод наименьших квадратов

В этом случае функциональный f (x), который будет решен,

$f (x) = {‖ C x + d ‖}_{2}^{2},$

возможно подвергните линейным ограничениям. Алгоритм генерирует строго выполнимый, выполняет итерации схождения, в пределе, к локальному решению. Каждая итерация включает приближенное решение большой линейной системы (порядка n, где n является длиной x). Матрицы итерации имеют структуру матричного C. В частности, метод предобусловленных методов сопряженных градиентов используется, чтобы приблизительно решить нормальные уравнения, i.e.,

$C^{T} C x = - C^{T} d,$

несмотря на то, что нормальные уравнения явным образом не формируются.

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

Функция умножения якобиана. lsqlin может решить задачу линейно-метода-наименьших-квадратов-с-ограничениями, не используя матричный C явным образом. Вместо этого это использует функцию умножения якобиана jmfun,

W = jmfun(Jinfo,Y,flag)

то, что вы обеспечиваете. Функция должна вычислить следующие продукты для матричного Y:

Если flag == 0 затем W = C'*(C*Y).
Если flag > 0 затем W = C*Y.
Если flag < 0 затем W = C'*Y.

Это может быть полезно, если C является большим, но содержит достаточно структуры, что можно записать jmfun не формируя C явным образом. Для примера смотрите функцию умножения якобиана с Линейным методом наименьших квадратов.

Метод Levenberg-Marquardt

Задача наименьших квадратов минимизирует функциональный f (x), который является суммой квадратов.

\min_{x} f (x) = {‖ F (x) ‖}_{2}^{2} = \sum_{i} F_{i}^{2} (x) .

(7)

Проблемы этого типа происходят в большом количестве практических применений, особенно те, которые включают подходящие функции модели к данным, такие как нелинейная оценка параметра. Этот проблемный тип также появляется в системах управления, где цель для выхода y (x,t), чтобы следовать за непрерывной траекторией модели φ (t) для векторного x и скалярного t. Эта проблема может быть описана как

\min_{x \in ℜ^{n}} \int_{t_{1}}^{t_{2}} {(y (x, t) - φ (t))}^{2} d t,

(8)

где y (x,t) и φ (t) является скалярными функциями.

Дискретизируйте интеграл, чтобы получить приближение

\min_{x \in ℜ^{n}} f (x) = \sum_{i = 1}^{m} {(y (x, t_{i}) - φ (t_{i}))}^{2},

(9)

где _ti равномерно расположен с интервалами. В этой проблеме векторный F (x)

$F (x) = [\begin{matrix} y (x, t_{1}) - φ (t_{1}) \\ y (x, t_{2}) - φ (t_{2}) \\ ... \\ y (x, t_{m}) - φ (t_{m}) \end{matrix}] .$

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

Обозначая m-by-n якобиевская матрица F (x) как J (x), вектор градиента из f (x) как G (x), матрица Гессиана f (x) как H (x) и матрица Гессиана каждого _Fi (x) как _Di (x),

\begin{array}{l} G (x) = 2 J (^{x) T} F (x) \\ H (x) = 2 J (^{x) T} J (x) + 2 Q (x), \end{array}

(10)

где

$Q (x) = \sum_{i = 1}^{m} F_{i} (x) \cdot D_{i} (x) .$

Свойство матричного Q (x) состоит в том, что, когда невязка ∥F (x) ∥ стремится по направлению к нулю, когда _xk приближается к решению, затем Q (x) также стремится по направлению к нулю. Так, когда ∥F (x) ∥ мал в решении, эффективный метод состоит в том, чтобы использовать направление Ньютона Гаусса в качестве базиса для процедуры оптимизации.

В каждой главной итерации k метод Ньютона Гаусса получает поисковое направление _dk, который является решением проблемы линейного метода наименьших квадратов

\min_{d_{k} \in ℜ^{n}} {‖ J (x_{k}) d_{k} + F (x_{k}) ‖}_{2}^{2} .

(11)

Направление, выведенное из этого метода, эквивалентно направлению Ньютона когда условия Q (x) = 0. Алгоритм может использовать поисковое направление _dk в качестве части стратегии поиска линии гарантировать, что функциональный f (x) уменьшается в каждой итерации.

Метод Ньютона Гаусса часто сталкивается с проблемами, когда термин второго порядка Q (x) ненезначителен. Метод Levenberg-Marquardt преодолевает эту проблему.

Метод Levenberg-Marquardt (см. [25] и [27]) использует поисковое направление, которое является решением линейной системы уравнений

(J {(x_{k})}^{T} J (x_{k}) + λ_{k} I) d_{k} = - J {(x_{k})}^{T} F (x_{k}),

(12)

или, опционально, уравнений

(J {(x_{k})}^{T} J (x_{k}) + λ_{k} diag (J {(x_{k})}^{T} J (x_{k}))) d_{k} = - J {(x_{k})}^{T} F (x_{k}),

(13)

где скалярный _λk управляет и величиной и направлением _dk и diag(A) означает матрицу диагональных терминов в A. Установите опцию ScaleProblem к 'none' выбрать Equation 12 или установить ScaleProblem к 'Jacobian' выбрать Equation 13.

Вы устанавливаете начальное значение параметра λ ₀ использований InitDamping опция. Иногда, 0.01 значение по умолчанию этой опции может быть неподходящим. Если вы находите, что алгоритм Levenberg-Marquardt делает мало начальных успехов, попробуйте установку InitDamping к различному значению от значения по умолчанию, такого как 1e2.

Когда _λk является нулем, направление_{, dk} идентичен тому из метода Ньютона Гаусса. Когда _λk стремится к бесконечности, _dk стремится к направлению наискорейшего спуска с величиной, стремящейся по направлению к нулю. Следовательно, для некоторого достаточно большого _λk, термин F (_xk + _dk) < F (_xk) сохраняется. Поэтому можно управлять термином _λk, чтобы гарантировать спуск, даже когда алгоритм сталкивается с условиями второго порядка, которые ограничивают КПД метода Ньютона Гаусса. То, когда шаг успешен (дает более низкое значение функции), алгоритм устанавливает λ _{k +1} = _λk/10. Когда шаг неудачен, алгоритм устанавливает λ _{k +1} = _λk *10.

Внутренне, алгоритм Levenberg-Marquardt использует допуск оптимальности (останавливающий критерий) 1e-4 времена функциональный допуск.

Метод Levenberg-Marquardt, поэтому, использует поисковое направление, которое является пересечением направления Ньютона Гаусса и направления наискорейшего спуска.

Другое преимущество для метода Levenberg-Marquardt состоит в том, когда якобиевский J имеет неполный ранг. В этом случае метод Ньютона Гаусса может иметь числовые проблемы, потому что проблема минимизации в уравнении 11 плохо изложена. В отличие от этого метод Levenberg-Marquardt имеет полный ранг в каждой итерации, и, поэтому, избегает этих проблем.

Следующий рисунок показывает итерации метода Levenberg-Marquardt при минимизации функции Розенброка, известно трудная проблема минимизации, которая находится в форме наименьших квадратов.

Метод Levenberg-Marquardt на функции Розенброка

Для большего количества полного описания этого рисунка, включая скрипты, которые генерируют итеративные точки, смотрите Банановую Минимизацию Функции.

Связанные ограничения в методе Levenberg-Marquardt

Когда проблема содержит связанные ограничения, lsqcurvefit и lsqnonlin измените итерации Levenberg-Marquardt. Если предложенная итеративная точка, x находится за пределами границ, алгоритм, проектирует шаг на самую близкую допустимую точку. Другими словами, с P, заданным как оператор projection, что проекты неосуществимые точки на выполнимую область, алгоритм изменяет предложенную точку x к P (x). По определению оператор проекции P работает с каждым _xi компонента независимо согласно

$P (x_{i}) = {\begin{array}{l} {lb}_{i} & если x_{i} < {lb}_{i} \\ {ub}_{i} & если x_{i} > {ub}_{i} \\ x_{i} & в противном случае \end{array}$

или, эквивалентно,

$P (x_{i}) = \min (\max (x_{i}, {lb}_{i}), {ub}_{i}) .$

Алгоритм изменяет останавливающийся критерий меры оптимальности первого порядка. Позвольте x быть предложенной итеративной точкой. В неограниченном случае останавливающийся критерий

{‖ \nabla f (x) ‖}_{\infty} \leq tol,

(14)

где tol является значением допуска оптимальности, которое является 1e-4*FunctionTolerance. В ограниченном случае останавливающийся критерий

{‖ x - P (x - \nabla f (x)) ‖}_{\infty}^{2} \leq tol {‖ \nabla f (x) ‖}_{\infty} .

(15)

Чтобы изучить этот критерий, сначала обратите внимание что, если x находится во внутренней части выполнимой области, то оператор P не оказывает влияния. Так, останавливающийся критерий становится

${‖ x - P (x - \nabla f (x)) ‖}_{\infty}^{2} = {‖ \nabla f (x) ‖}_{\infty}^{2} \leq tol {‖ \nabla f (x) ‖}_{\infty},$

который совпадает с исходным неограниченным критерием остановки, ${‖ \nabla f (x) ‖}_{\infty} \leq tol$ . Если граничное ограничение активно, означая, что x – ∇f (x) не выполним, то в точке, где алгоритм должен остановиться, градиент в точке на контуре перпендикулярен контуру. Поэтому точка x равна P (x – ∇f (x)), проекция шага наискорейшего спуска, как показано в следующем рисунке.

Levenberg-Marquardt останавливающееся условие

Sketch of x minus the projection of x minus gradient of f(x)

Документация