Оптимальная аппроксимация сплайнами

в 17:42, , рубрики: data mining, python3, Алгоритмы, математика, машинное обучение, оптимизация, сплайны, метки:

Пусть нам дан набор точек Оптимальная аппроксимация сплайнами - 1 и соответствующий им набор положительных весов Оптимальная аппроксимация сплайнами - 2. Мы считаем, что некоторые точки могут быть важнее других (если нет, то все веса одинаковые). Неформально говоря, мы хотим, чтобы на соответствующем интервале была проведена красивая кривая таким образом, чтобы она «лучше всего» проходила через эти данные.

Оптимальная аппроксимация сплайнами - 3

Под катом находится алгоритм, раскрывающий, каким образом сплайны позволяют строить подобную красивую регрессию:
Оптимальная аппроксимация сплайнами - 4

Основные определения

Функция s(x) на интервале [a, b] называется сплайном степени k на сетке с горизонтальными узлами Оптимальная аппроксимация сплайнами - 5, если выполняются следующие свойства:

  • На интервалах Оптимальная аппроксимация сплайнами - 6 функция s(x) является полиномом k-й степени.
  • n-ая производная функции s(x) непрерывна в любой точке [a, b] для любого n = 1,…, k-1.

Заметим, что для построения сплайна нужно для начала задать сетку из горизонтальных узлов. Расположим их таким образом, чтобы внутри интервала (a, b) стояло g узлов, а по краям — k+1: Оптимальная аппроксимация сплайнами - 7 и Оптимальная аппроксимация сплайнами - 8.

Каждый сплайн в точке Оптимальная аппроксимация сплайнами - 9 может быть представлен в базисной форме:

Оптимальная аппроксимация сплайнами - 10

где Оптимальная аппроксимация сплайнами - 11B-сплайн k+1-го порядка:

Оптимальная аппроксимация сплайнами - 12

Оптимальная аппроксимация сплайнами - 13

Оптимальная аппроксимация сплайнами - 14

Вот как, например, выглядит базис на сетке из g = 9 узлов, равномерно распределенных на интервале [0, 1]:
Оптимальная аппроксимация сплайнами - 15
Сходу разобраться в построении сплайнов через B-сплайны очень сложно. Больше информации можно найти здесь: www.brnt.eu/phd/node11.

Аппроксимация с заданными горизонтальными узлами

Итак, мы выяснили что сплайн определяется однозначно узлами и коэффициентами. Допустим, что узлы Оптимальная аппроксимация сплайнами - 16 нам известны. Также на вход подается набор данных Оптимальная аппроксимация сплайнами - 17 с соответствующими весами Оптимальная аппроксимация сплайнами - 18. Необходимо найти коэффициенты Оптимальная аппроксимация сплайнами - 19, максимально приближающие кривую сплайна к данным. Строго говоря, они должны доставлять минимум функции

Оптимальная аппроксимация сплайнами - 20

Для удобства запишем в матричном виде:

Оптимальная аппроксимация сплайнами - 21

где

Оптимальная аппроксимация сплайнами - 22

Оптимальная аппроксимация сплайнами - 23

Заметим, что матрица E блочно-диагональная. Минимум достигается когда градиент ошибки по коэффициентам будет равен нулю:

Оптимальная аппроксимация сплайнами - 24

Зададим оператор Оптимальная аппроксимация сплайнами - 25, обозначающий взвешенное скалярное произведение:

Оптимальная аппроксимация сплайнами - 26

Оптимальная аппроксимация сплайнами - 27

Пусть также

Оптимальная аппроксимация сплайнами - 28

Оптимальная аппроксимация сплайнами - 29

Тогда вся задача и все предыдущие формулы сводятся к решению простой системы линейных уравнений:

Оптимальная аппроксимация сплайнами - 30

где матрица А (2k+1)-диагональная, так как Оптимальная аппроксимация сплайнами - 31, если |i — j| > k. Также матрица А симметричная и положительно-определенная, следовательно решение возможно быстро найти с помощью разложения Холецкого (существует также алгоритм для разреженных матриц).
И вот, решая систему, получаем желаемый результат:
Оптимальная аппроксимация сплайнами - 32

Сглаживание

Однако, далеко не всегда все так хорошо. При малом количестве данных по отношению к количеству узлов и степени сплайна может возникнуть проблема т.н. сверхподгонки (overfitting). Вот пример «плохого» кубического сплайна, при этом идеально проходящего сквозь данные:
Оптимальная аппроксимация сплайнами - 33
Окей, кривая уже не такая уж и красивая. Попытаемся уменьшить так называемые колебания сплайна. Для этого мы попробуем «сгладить» его k-ю производную. Другими словами, мы минимизируем разницу между производной слева и производной справа от каждого узла:

Оптимальная аппроксимация сплайнами - 34

Разложив сплайн в базисную форму, мы получаем:

Оптимальная аппроксимация сплайнами - 35

Оптимальная аппроксимация сплайнами - 36

Давайте рассмотрим ошибку

Оптимальная аппроксимация сплайнами - 37

Здесь q — вес функции, влияющей на сглаживание, и

Оптимальная аппроксимация сплайнами - 38

Оптимальная аппроксимация сплайнами - 39

Новая система уравнений:

Оптимальная аппроксимация сплайнами - 40

где

Оптимальная аппроксимация сплайнами - 41

Ранг матрицы B равен g. Она симметричная и, так как q > 0, A + qB будет положительно определенной. Поэтому разложение Холецкого по-прежнему применимо к новой системе уравнений. Однако, матрица B вырожденная и при слишком больших значениях q могут возникнуть численные ошибки.
При совсем маленьком значении q = 1e-9 вид кривой изменяется очень слабо.
Оптимальная аппроксимация сплайнами - 42
Но при q = 1e-7 в данном примере уже достигается достаточное сглаживание.
Оптимальная аппроксимация сплайнами - 43

Аппроксимация с неизвестными горизонтальными узлами

Представим теперь, что задача такая же, как и прежде, за исключением того, что мы не знаем как узлы расположены на сетке. На вход кроме данных подается только количество узлов g, интервал [a, b] и степень сплайна k. Попробуем наивно предположить, что лучше всего расположить узлы равномерно на интервале:
Оптимальная аппроксимация сплайнами - 44

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

Оптимальная аппроксимация сплайнами - 45

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

Оптимальная аппроксимация сплайнами - 46

Положительный параметр p — вес штрафной функции. Чем больше его значение, тем быстрее узлы будут удаляться друг от друга и стремиться к равномерному расположению.
Для решения данной задачи мы используем метод сопряженных градиентов. Его прелесть заключается в том, что для квадратичной функции он сходится за фиксированное (в данном случае g) количество шагов.

  1. Инициализируем направление Оптимальная аппроксимация сплайнами - 47.
    Как рассчитать производную ошибки по узлам?

    Производная суммы квадратов по узлу:

    Оптимальная аппроксимация сплайнами - 48

    Для того, чтобы рассчитать влияние положения узла на значения сплайна, нужно рассмотреть B-сплайны Оптимальная аппроксимация сплайнами - 49 на новых узлах Оптимальная аппроксимация сплайнами - 50 и с новыми коэффициентами Оптимальная аппроксимация сплайнами - 51

    Оптимальная аппроксимация сплайнами - 52

    Оптимальная аппроксимация сплайнами - 53

    Оптимальная аппроксимация сплайнами - 54

    Производная штрафной функции:

    Оптимальная аппроксимация сплайнами - 55

    На производную функции сглаживания без слез не взглянешь:

    Оптимальная аппроксимация сплайнами - 56
  2. Для j = 0,…, g-1
    • Зададим функцию
      Оптимальная аппроксимация сплайнами - 57

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

    • Обновляем значения узлов:
      Оптимальная аппроксимация сплайнами - 58
    • Обновляем вектор направления:
      Оптимальная аппроксимация сплайнами - 59
  3. Если
    Оптимальная аппроксимация сплайнами - 60

    и

    Оптимальная аппроксимация сплайнами - 61

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

Решение задачи одномерной минимизации

Для того, чтобы найти значение Оптимальная аппроксимация сплайнами - 62, доставляющее минимум функции

Оптимальная аппроксимация сплайнами - 63

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

  1. Пусть первая и последняя компоненты вектора направления равны нулю: Оптимальная аппроксимация сплайнами - 65. Зададим также максимально возможный шаг вдоль этого направления:
    Оптимальная аппроксимация сплайнами - 66

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

  2. Инициируем k = 0 и начальные шаги: Оптимальная аппроксимация сплайнами - 67, Оптимальная аппроксимация сплайнами - 68, Оптимальная аппроксимация сплайнами - 69.
  3. До тех пор, пока Оптимальная аппроксимация сплайнами - 70:
    • Задаём
      Оптимальная аппроксимация сплайнами - 71

      и уменьшаем шаг

      Оптимальная аппроксимация сплайнами - 72
    • k = k + 1

  4. Если k > 0, то возвращаем α* = α1. Иначе:
    • До тех пор, пока Оптимальная аппроксимация сплайнами - 73: α0 = α1, α1 = α2 и
      Оптимальная аппроксимация сплайнами - 74

    • Возвращаем Оптимальная аппроксимация сплайнами - 75, где Оптимальная аппроксимация сплайнами - 76 — корень уравнения I'(α) = 0 и I(α) — аппроксимация функции ошибки:
      Оптимальная аппроксимация сплайнами - 77

      где

      Оптимальная аппроксимация сплайнами - 78

      Оптимальная аппроксимация сплайнами - 79

Коэффициенты ai и bi могут быть найдены из уравнений

Оптимальная аппроксимация сплайнами - 80

и

Оптимальная аппроксимация сплайнами - 81

Объяснение алгоритма:
Идея заключается в том, чтобы расставить три точки α0 < α1 < α2 таким образом, чтобы по значениям ошибок, достигаемых в этих точках, можно было построить простую аппроксимирующую функцию и вернуть её минимум. Притом значение ошибки в α1 должно быть меньше, чем значение ошибки в α0 и α2.

Находим начальное приближение α1 из условия S'(α1)=0, где S(α) — функция вида

Оптимальная аппроксимация сплайнами - 82

Константы c0 и c1 находятся из условий Оптимальная аппроксимация сплайнами - 83 и Оптимальная аппроксимация сплайнами - 84.

Если мы просчитались с начальным приближением, то мы уменьшаем шаг α1 до тех пор, пока он доставляет большее значение ошибки, чем α0. Выбор Оптимальная аппроксимация сплайнами - 85 исходит из условия Оптимальная аппроксимация сплайнами - 86, где Q(α) — парабола интерполирующая функцию ошибки Оптимальная аппроксимация сплайнами - 87: Оптимальная аппроксимация сплайнами - 88, Оптимальная аппроксимация сплайнами - 89 и Оптимальная аппроксимация сплайнами - 90.
Если k > 0, то мы нашли значение α1, такое что при его выборе значение ошибки будет меньше, чем при выборе α0 и α2, и мы возвращаем его в качестве грубого приближения α*.

Если же наше первоначальное приближение было верным, то мы пытаемся найти шаг α2, такой что Оптимальная аппроксимация сплайнами - 91. Он будет найден между α1 и αmax, так как αmax — точка сингулярности для штрафной функции.

Когда найдены все три значения α0, α1 и α2, мы представляем функцию ошибки в виде суммы двух функций, приближающих разность квадратов и функцию штрафа. Функция Q(α) — парабола, чьи коэффициенты могут быть найдены, так как мы знаем её значения в трех точках. Функция R(α) уходит на бесконечность при α, стремящемся к αmax. Коэффициенты bi также могут быть найдены из системы из трех уравнений. В результате, мы приходим к уравнению, которое может быть приведено к квадратному и легко решено:

Оптимальная аппроксимация сплайнами - 92

И вот, для сравнения, результат оптимально построенного сплайна:
Оптимальная аппроксимация сплайнами - 93

Ну и для тех, кому может пригодиться: реализация на Python.

Автор: The_Freeman

Источник

* - обязательные к заполнению поля


https://ajax.googleapis.com/ajax/libs/jquery/3.4.1/jquery.min.js