Главная       Скачать       Коммерческая поддержка       FAQ       Forum       О нас       Английская версия

Аппроксимация линейным или нелинейным МНК

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

Содержание

    1 Виды аппроксимации
           Линейный МНК и родственные задачи
           Нелинейный МНК
           Рекомендации по выбору
    2 Аппроксимация полиномами
           Получение барицентрического представления
           Конверсия в степенной базис
           Примеры
    3 Аппроксимация рациональными функциями
           Получение барицентрического представления
    4 Аппроксимация кубическими сплайнами
           Обзор алгоритмов
           Регрессионный сплайн со штрафной функцией
           Достоинства
           Настройка
           Практический пример
    5 Аппроксимация с использованием линейного МНК
           Задачи без ограничений
           Задачи с линейными ограничениями
           Примеры
    6 Аппроксимация функциями общего вида
           Общие принципы
           Обзор методов
           Good fit vs .bad fit
           Выбор критериев остановки
           Масштабируя переменные
           Добавляя ограничения
           Дополнительные настройки
           Примеры

Виды аппроксимации

Линейный МНК и родственные задачи

Большинство методов аппроксимации, входящих в состав ALGLIB, в конечном итоге сводится к решению проблемы линейного МНК:

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

Нелинейный МНК

Пакет ALGLIB поддерживает нелинейную аппроксимацию с использованием определенной пользователем функции. Для поиска оптимума используется метод Левенберга-Марквардта, входящий в состав пакета ALGLIB. Интерфейс, рассматриваемый на этой странице, фактически является дружественной к пользователю оберткой вокруг лежащего в основе метода оптимизации.

Рекомендации по выбору

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

Аппроксимация полиномами

Получение барицентрического представления

Полиномиальная аппроксимация осуществляется с использованием подпрограмм polynomialfit (аппроксимация без индивидуальных ограничений и весовых коэфициентов) и polynomialfitwc (аппроксимация с ограничениями и/или индивидуальными весовыми коэффициентами). Поддерживается произвольное количество ограничений на значение функции - f(xc)=yc - или на значение первой производной - df(xc)/dx=yc. Трудоемкость решения составляет O(N·M 2) (здесь N - количество точек, M - размер базиса).

Замечание #1
Избыточные ограничения могут быть несовместны, особенно в случае аппроксимации полиномом низкой степени. Дополнительная информация по данному вопросу: см. 'Reference manual', комментарии к подпрограмме polynomialfitwc.

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

Конверсия в степенной базис

Пакет ALGLIB использует барицентрическое представление для операций с полиномами, потому что оно характеризуется большей стабильностью, чем представление в степенном базисе (в виде суммы степеней x). Результатом работы функций является объект barycentricinterpolant, хранящий коэффициенты барицентрической модели. Субпакеты polint и ratint содержат ряд функций для работы с этим объектом. Дополнительная информация доступна на странице User Guide, посвященной полиномальной интерполяции.

Однако две функции наиболее часто оказываются востребованными при полиномиальной аппроксимации. Иногда требуется получить именно коэффициенты полинома, а не просто "черный ящик" (барицентрический интерполянт), позволяющий вычислить значение полинома в любой точке. Для этого можно использовать следующие функции, позволяющие конвертировать барицентрический интерполянт в другую форму:

Примеры

Следующие примеры посвящены использованию полиномов для решения задач аппроксимации. Для ознакомления с полиномиальной интерполяцией вообще мы рекомендуем ознакомиться с описанием субпакета polint в Reference Manual или прочитать соответствующую страницу User Guide.

Аппроксимация рациональными функциями

Получение барицентрического представления

Пакет ALGLIB поддерживает рациональную аппроксимацию с использованием базиса Флоатера-Хорманна. Аппроксимация осуществляется двумя подпрограммами. Подпрограмма barycentricfitfloaterhormann решает задачу аппроксимации без ограничений и индивидуальных весов. Результатом работы является структура barycentricinterpolant, хранящая барицентрическое представление рациональной функции. Пакет ALGLIB содержит ряд функций для работы с этим объектом, описанных в документации к субпакету ratint и на странице User Guide, посвященной рациональной интерполяции.

Подпрограмма barycentricfitfloaterhormannwc решает более сложную задачу: аппроксимацию с индивидуальными весовыми коэффициентами и ограничениями на значения построенной модели. Поддерживается произвольное количество ограничений вида f(xc)=yc или df(xc)/dx=yc.

Замечание #2
Следует учитывать, что избыточные ограничения могут оказаться несовместимы. Более подробно этот вопрос рассмотрен в комментариях к подпрограмме barycentricfitfloaterhormannwc.

Обе подпрограммы используют базис Флоатера-Хорманна, построенный на равномерной сетке с M узлами. ункции перебирают разные значения D в интервале [0,9] в поисках оптимального, т.е. дающего минимальную взвешенную ошибку. Поскольку в такой формулировке задача линейна, то для решения используется линейный МНК-солвер и трудоемкость решения равна O(N*M^2), где N - число точек, M - размерность базиса.

Важно понимать, что для аппроксимации используется множество рациональных функций, гарантированно не имеющих полюсов - т.е. подмножество множества рациональных функций. С одной стороны, это делает невозможными резкие осцилляции интерполянта. С другой - снижает точность полученных моделей. Базис Флоатера-Хорманна позволяет аппроксимировать гладкую функцию с точностью 10 -3...10 -5, используя умеренное количество точек, однако для дальнейшего увеличения точности может потребоваться очень большое количество точек.

Аппроксимация кубическими сплайнами

Обзор алгоритмов

Пакет ALGLIB поддерживает три типа аппроксимации сплайнами:

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

Регрессионный сплайн со штрафной функцией

Регрессионный сплайн со штрафной функцией - это современный алгоритм аппроксимации, пригодный для работы с зашумленными экспериментальными данными. Это сплайн с M равномерно распределенными узлами, коэффициенты которого получаются путем минимизации ошибки аппроксимации LS и штрафной функции P, подавляющей нелинейность:

Для построения регрессионного сплайна можно использовать две функции:

Обе функции сводят задачу минимизации LS+P к линейной проблеме МНК. Трудоемкость решения составляет O(N·M 2+M 3). Результатом является объект spline1dinterpolant, хранящая кубический сплайн в упакованном виде. Пакет ALGLIB содержит ряд функций для работы с этим объектом, описанных в документации к субпакету spline1d и на странице User Guide, посвященной интерполяции кубическими сплайнами.

Достоинства

Ключевыми чертами регрессионного сплайна являются:

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

Настройка

Как и все алгоритмы, регрессионный сплайн нуждается в настройке, однако она осуществляется очень легко и просто:

  1. выберите количество узлов (степеней свободы) M. Количество узлов должно быть достаточно большим, чтобы обеспечить требуемую гибкость, и достаточно малым, чтобы обеспечить требуемое быстродействие. Не бойтесь выбрать слишком большое значение - штрафная функция помешает сплайну стать "слишком гибким". Можете сразу начать с сотни узлов, и увеличивать это число на сотню, до тех пор, пока не решите, что достаточно. Если быстродействие не является проблемой - можете попробовать сделать M=3·N.
  2. выберите коэффициент регуляризации ρ. Чем больше ρ, тем больше степень сглаживания. При установке максимальной степени сглаживания вы получите прямую линию (не обязательно горизонтальную). Этот коэффициент автоматически масштабируется таким образом, что рабочий диапазон его значений равен [-15,+15], а в большинстве практических задач оптимальное ρ находится в диапазоне [-2,+6].

Таким образом, настройка алгоритма состоит из двух независимых стадий, практически не связаных друг с другом. Каждый из параметров настраивается отдельно от другого, при этом основное внимание должно уделаться настройке ρ, а количество степеней свободы M лучше настраивать по принципу "установи и забудь".

Практический пример

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

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

На втором графике показано три сплайна, соответствующих различным значениям ρ. Синяя кривая соответствует очень низкому уровню сглаживания (хотя не самому низкому из возможных - ρ может быть и отрицательным). Темно-красная кривая соответствует умеренному уровню сглаживания, достаточному для подавления высокочастотной шумовой составляющей. Наконец, зеленая кривая показывает нам результат, полученный с использованием ρ=6.0 - можно видеть, что шумовая составляющая была полностью подавлена.

Замечание #3
Обратите внимание, что число степеней свободы M намного больше числа точек N, но какие-либо признаки переобучения отсутствуют.

Замечание #4
При дальнейшем увеличении ρ наш сплайн стал бы постепенно приближаться к прямой линии.

Исходный код примера, показанного выше, недоступен для скачивания, потому что он несовместим с используемой нами системой автоматической генерации примеров. Однако ALGLIB reference manual содержит исходный код других примеров аппроксимации кубическими сплайнами:

Аппроксимация с использованием линейного МНК

Задачи без ограничений

Пакет ALGLIB содержит две функции для решения задач линейного МНК без ограничений.

Линейный МНК может быть применен, если аппроксимирующая функция представляется, как линейная комбинация базисных функций. При этом сами базисные функции могут быть нелинейны по x.

Например, первая из двух функций, приведенных ниже, представляет аппроксимирующий полином, как сумму степеней x. Вторая функция служит для решения другой проблемы - аппроксимации данных как суммы синусов.

Решением задачи является вектор c, минимизирующий E(x) - общую ошибку аппроксимации. Ниже приведено два варианта E(x) - для задачи с весовыми коэффициентами и без:

Замечание #5
Обратите внимание, что весовые коэффициенты, если они есть, возводятся в квадрат.

В такой формулировке задача сводится к решению системы линейных уравнений относительно коэффициентов c, что позволяет использовать линейные солверы, находящие ответ за фиксированное время порядка O(N·M 2). Таким образом, вам не придется ждать неопределенное время, пока итеративный алгоритм сойдется к ответу. Другим достоинством, вытекающим из линейности задачи, является то, что от вас не требуется многократно вычислять значения базисных функций в различных точках. Достаточно один раз вычислить матрицу значений базисных функций и передать её в алгоритм.

Замечание #6
Для решения используется основанный на SVD солвер, позволяющий решать плохо обусловленные проблемы и даже проблемы с вырожденной (линейно зависимой) системой базисных функций.

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

Задачи с линейными ограничениями

Другим важным видом задач МНК являются задачи, где на вектор параметров c накладываются линейные ограничения в виде равенств:

Для решения таких задач можно использовать две функции:

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

Замечание #8
Ограничения могут быть несовместными, т.е. может не существовать такого c, которое удовлетворяет им всем. В этом случае алгоритм вернет соответствующий код ошибки.

Примеры

Этот раздел содержит ссылки на примеры применения линейного МНК:

Аппроксимация функциями общего вида

Общие принципы

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

Нелинейная аппроксимация очень существенно отличается от линейной:

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

  1. сначала мы создаем объект-оптимизатор при помощи одной из функций-конструкторов
  2. мы настраиваем оптимизатор, задавая условия остановки (и меняя другие параметры, при необходимости)
  3. мы запускаем процесс решения, вызывая функцию lsfitfit
  4. после завершения процесса мы получаем результат при помощи функции lsfitresults

Обзор методов

Пользователи могут использовать три основных варианта алгоритма:

Каждый из алгоритмов доступен в двух вариантах - взвешенном, при котором каждой точке назначен индивидуальный весовой коэффициент, и невзвешенном, при котором все точки имеют равный вес. Буквы в названии варианта алгоритма являются суффиксом, который дописывается к имени подпрограммы lsfitcreate, использующейся для создания оптимизатора. Если используется взвешенный вариант, дополнительно приписывается суффикс W. Итого мы получаем шесть вариантов:

Какую же схему следует выбрать? Для быстрого старта мы рекомендуем начать со схемы F, потому что она наиболее проста в использовании. Вы просто пишете код, вычисляющий значение функции, а пакет ALGLIB берет на себя вопросы, связанные с численным дифференцированием.

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

Good fit vs .bad fit

Известно, что метод Левенберга-Марквардта быстрее всего сходится, если все точки лежат близко от кривой (good fit). Если задача относится к категории "bad fit", то скорость сходимости снижается в несколько раз. Если это является проблемой, вы можете попробовать использовать FGH-вариант алгоритма. Он требует меньше итераций, но стоимость каждой итерации выше, так что не всегда понятно, окупится он или нет. В тех редких случаях, когда вам доступен дешевый Гессиан, FGH-вариант скорее всего окупит себя.

Выбор критериев остановки

Пакет ALGLIB предлагает пользователям три критерия остановки:

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

Мы рекомендуем использовать первый критерий (достаточно малый шаг). Выбор величины шага зависит от решаемой задачи, однако мы рекомендуем брать очень маленькое значение, намного меньше, чем точность, с которой вы хотите найти решение. Это связано с тем, что иногда, в местах со сложным рельефом, алгоритм может делать очень маленькие шаги. Т.е если длина последнего шага равна 0.001, то это не значит, что расстояние до минимума имеет тот же порядок. Например, это может значить, что наша квадратичная модель устарела и алгоритм как раз пытается понять, не следует ли обновить её. Или что мы прошли сложный поворот по дну оврага и ещё не набрали скорость перед выходом на финишную прямую.

Замечание #9
Дополнительно вы можете установить ограничение по количеству шагов - если хотите гарантировать, что алгоритм не проработает слишком долго. Например, такое условие имеет смысл во внедренных системах, когда решение должно быть получено достаточно оперативно, либо не получено вообще.

Замечание #10
критерий "малое значение градиента" недоступен пользователю, т.к. сложно понять, какое же значение "градиента взвешенной суммы квадратов отклонений" будет "достаточно малым" для данной проблемы. Критерии, основанные на длине шага или снижении функции, более интуитивны и понятны.

Масштабируя переменные

Перед началом работы рекомендуется установить масштаб параметров при помощи функции lsfitsetscale. Без установки масштаба можно обойтись, если ваша проблема хорошо отмасштабирована. Если величина некоторых параметров различается более чем в сто раз, рекомендуется установить масштаб. При более значительном различии (от тысячи раз) установка масштаба может быть необходимой для корректной работы оптимизатора.

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

Добавляя ограничения

Пакет ALGLIB поддерживает аппроксимацию с простыми ограничениями на параметры, т.е. ограничениями вида l≤c≤u. Простые ограничения устанавливаются при помощи функции lsfitsetbc. Эти ограничения обрабатываются очень эффективно - N ограничений требует всего лишь O(N) дополнительных операций при каждом вычислении функции. Наконец, эти ограничения всегда выполняются в точности. Мы не вычисляем функцию в точках, координаты которых лежат за пределами интервала, задаваемого ограничениями. Результат оптимизации также всегда строго удовлетворяет ограничениям (находится либо во внутренней области, либо на границе).

Дополнительные настройки

В дополнение к указанным выше настройкам пользователь может:

Примеры

Этот раздел содержит ссылки на примеры применения нелинейного МНК:

This article is intended for personal use only.

Скачать ALGLIB

C#

Исходный код на C#

Downloads page

 

C++

Исходный код на C++

Downloads page

 

C++, арифметика высокой точности

Исходный код на C++, использующий библиотеки MPFR/GMP.

Исходный код GMP доступен на сайте gmplib.org. Исходный код MPFR доступен на сайте www.mpfr.org.

Downloads page

 

FreePascal

Исходный код на Free Pascal.

Downloads page

 

Delphi

Исходный код на Delphi.

Downloads page

 

VB.NET

Исходный код на VB.NET.

Downloads page

 

VBA

Исходный код на VBA.

Downloads page

 

Python

Исходный код на Python (CPython и IronPython).

Downloads page

 

 

ALGLIB® - numerical analysis library, 1999-2014.
ALGLIB is registered trademark of the ALGLIB Project.