Алгоритмы наименьших квадратов (Model Fitting)

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

Наименьшие квадраты, в целом, является проблемой нахождения вектора 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 (<reservedrangesplaceholder0>)	Ограничения
`mldivide`	C · x – d	Ничего
`lsqnonneg`	C · x – d	<reservedrangesplaceholder0> ≥ 0
`lsqlin`	C · x – d	Граница, линейная
`lsqnonlin`	Общие F (x)	Связанный
`lsqcurvefit`	F (x, xdata) – ydata	Связанный

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

lsqlin interior-point
lsqlin active-set
Trust- области отражающий (нелинейный или линейный методом наименьших квадратов)
Левенберг-Марквардт (нелинейные методы наименьших квадратов)
Алгоритм, используемый lsqnonneg

Все алгоритмы, кроме lsqlin активный набор являются масштабными; см. «Алгоритмы большого и среднего масштаба». Общее исследование нелинейных методов наименьших квадратов смотрите в Dennis [8]. Конкретные детали метода Левенберга-Марквардта можно найти в Moré [28].

Линейные методы наименьших квадратов: Interior-Point или Active-Set

lsqlin 'interior-point' алгоритм использует Алгоритм внутренней точки выпуклого квадратиков, и 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 равной 2C^TC и вектор c to (-2C^Tг). (Аддитивный термин d^Td не влияет на положение минимума.) После этого переформулирования lsqlin задача, quadprog вычисляет решение.

Примечание

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

Trust- Области светоотражающий метод наименьших квадратов

Алгоритм Trust- Области Reflective методом наименьших квадратов

Многие методы, используемые в решателях 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 1 и _s 2, где s 1 находится в направлении градиента g, и s 2 является либо приблизительным направлением Ньютона, т.е. решением

H \cdot s_{2} = - g,

(3)

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

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

(4)

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

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

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

Эти четыре шага повторяются до сходимости. Размерность Β доверительной области регулируется согласно стандартным правилам. В частности, это уменьшено, если шаг испытания не принят, т.е. f (x + s) ≥ <reservedrangesplaceholder1> (<reservedrangesplaceholder0>). Для обсуждения этого аспекта см. [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) не используются.

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

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

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

Система больших Шкал Линейного метода наименьших квадратов

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

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

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

$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 явным образом. Для получения примера см. Функцию умножения якобиана с Линейным методом наименьших квадратов.

Метод Левенберга-Марквардта

Задача наименьших квадратов минимизирует 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 n якобиевской матрицей F (<reservedrangesplaceholder14>) как J (<reservedrangesplaceholder12>), вектор градиента f _{(<reservedrangesplaceholder10>)} как G (<reservedrangesplaceholder8>), матрицей Мешковины f (<reservedrangesplaceholder6>) как H (<reservedrangesplaceholder4>) и матрицей Мешковины каждого Fi (<reservedrangesplaceholder2>) как Di (<reservedrangesplaceholder0>),

\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) неотрицателен. Метод Левенберга-Марквардта преодолевает эту проблему.

Метод Левенберга-Марквардта (см. [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' выбрать уравнение 12 или задать ScaleProblem на 'Jacobian' для выбора уравнения 13.

Начальное значение параметра вы устанавливаете λ 0 используя InitDamping опция. Иногда, 0.01 значение по умолчанию этой опции может быть непригодным. Если вы обнаружили, что алгоритм Левенберга-Марквардта мало прогрессирует, попробуйте задать InitDamping к значению, отличному от значения по умолчанию, такому как 1e2.

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

Внутренне алгоритм Левенберга-Марквардта использует допуск оптимальности (критерий остановки) 1e-4 умножает функцию на допуск.

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

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

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

Метод Левенберга-Марквардта о функции Розенбрка

Более полное описание этого рисунка, включая скрипты, которые генерируют итерационные точки, смотрите в Banana Function Minimization.

Связанные ограничения в методе Левенберга-Марквардта

Когда задача содержит связанные ограничения, lsqcurvefit и lsqnonlin измените итерации Левенберга - Марквардта. Если предложенная итерационная 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)), проекции шага наискорейшего спуска, как показано на следующем рисунке.

Левенберг-Марквардт Остановка Условия

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

См. также

lsqcurvefit | lsqlin | lsqnonlin | lsqnonneg | quadprog

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