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

С математической точки зрения кривая, заданная вершинами многоугольника, зависит от интерполяции или аппроксимации, устанавливающей связь кривой и многоугольника. Здесь основой является выбор базисных функций. Как было отмечено в разд. 5-8, базис Бернштейна порождает кривые Безье вида (5-62), но он обладает двумя свойствами, которые ограничивают гибкость кривых. Во-первых, количество вершин многоугольника жестко задает порядок многочлена. Например, кубическая кривая должна быть задана четырьмя вершинами и тремя отрезками. Многоугольник из шести точек всегда порождает кривую пятого порядка. Единственный способ понизить степень кривой - это сократить количество вершин, а повысить степень кривой - увеличить их число.

Второе ограничение следует из глобальной природы базиса Бернштейна. Это означает, что величина аппроксимирующих функций из уравнения (5-63) ненулевая для всех значений параметра на кривой. Любая точка на кривой Безье зависит от всех определяющих вершин, поэтому изменение какой-либо одной вершины оказывает влияние на всю кривую. Локальные воздействия на кривую невозможны.

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

Существует неглобальный базис, называемый базисом В-сплайна, включающий базис Бернштейна как частный случай. В-сплайны неглобальны, так как с каждой вершиной связана своя базисная функция. Поэтому влияние каждой вершины на кривую проявляется только при тех значениях параметра, где соответствующая базисная функция не равна нулю. Базис В-сплайна также позволяет менять порядок базисных функций и, следовательно, всей кривой без изменения количества вершин. Теория В-сплайнов была предложена в работе . Рекурсивное определение для численного решения было выведено независимо Коксом и де Буром . Гордон и Ризенфельд , определяли кривые через базис В-сплайна.

Пусть определяет кривую как функцию от параметра , тогда В-сплайн имеет вид

, , , (5-83)

где есть вершина многоугольника, a - нормализованные функции базиса В-сплайна.

Для -й нормализованной функции базиса порядка (степени ) функции базиса определяются рекурсивными формулами Кокса-де Бура:

(5-84а)

Величины - это элементы узлового вектора, удовлетворяющие отношению . Параметр изменяется от до вдоль кривой . Считается, что .

Формально В-сплайн определяется как полиномиальный сплайн порядка (степени ), так как он удовлетворяет следующим условиям:

Функция является полиномом степени на каждом интервале .

И ее производные порядка непрерывны вдоль всей кривой.

Так, например, В-сплайн четвертого порядка - это кусочная кубическая кривая.

Из того что В-сплайн задается базисом В-сплайна, сразу следует еще несколько его свойств:

Сумма базисных функций В-сплайна для любого значения параметра (см. работы и )

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

Кроме , все базисные функции имеют ровно один максимум.

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

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

Общая форма кривой повторяет форму определяющего многоугольника.

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

Кривая лежит внутри выпуклой оболочки определяющего многоугольника.

Последнее свойство В-сплайна сильнее, чем у кривых Безье. У В-сплайна порядка (степени ) точки кривой лежат внутри выпуклой оболочки соседних точек. Таким образом, все точки на В-сплайне должны лежать внутри объединения всех выпуклых оболочек последовательных вершин. На рис. 5-32 приводится иллюстрация для различных значений , причем выпуклые оболочки выделены серым цветом. В частности, при выпуклая оболочка совпадает с многоугольником, т. е. В-сплайн - это сам многоугольник.

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

Если имеется совпадающих вершин, т.е. , то выпуклая оболочка вершин от до - это сама вершина. Отсюда следует, что В-сплайн должен проходить через вершину . На рис. 5-34 изображен пример такой точки для . Далее, так как В-сплайн везде непрерывен, он также непрерывен в .

Наконец, заметим, что свойство непрерывности плавно переводит В-сплайн во вложенные отрезки прямой, как показано на рис. 5-35.

Уравнения (5-84) указывают, что выбор узлового вектора оказывает существенное влияние на базисные функции В-сплайна и, следовательно, на сам В-сплайн. Единственное требование к узловому вектору: , т. е. это монотонно возрастающая последовательность вещественных чисел. Обычно используются три типа узловых векторов: равномерные, открытые равномерные (или открытые) и неравномерные.

Отдельные узловые значения равномерного узлового вектора распределены на одинаковом расстоянии, например

В частности, равномерные узловые векторы обычно начинаются в нуле и увеличиваются на 1 к некоторому максимальному значению или нормируются в диапазоне между 0 и 1 равными десятичными значениями, например,

Рис. 5-32 Свойства выпуклой оболочки В-сплайнов.

Рис. 5-33 Свойства выпуклой оболочки В-сплайнов для коллинеарных сегментов кривой. (а) Внутренние вершины определяющего многоугольника; (b) вершины в конце определяющего многоугольника.

Для данного порядка равномерные узловые векторы порождают периодические равномерные функции базиса, для которых

То есть каждая функция базиса - это параллельный перенос другой функции, см. рис. 5-36.

У открытого равномерного узлового вектора количество одинаковых узловых значений в концах равно порядку базисной функции В-сплайна. Внутренние узловые значения распределены равномерно.

Рис. 5-34 Выпуклая оболочка для кратных вершин, .

Рис. 5-35 Плавное превращение в отрезки прямой.

Рис. 5-36 Базисные функции периодического равномерного В-сплайна, , , .

Несколько примеров с целыми приращениями:

,

или для нормализованных приращений

Рис. 5-37 Базисные функции открытого равномерного В-сплайна, , , .

Формально открытый равномерный узловой вектор определяется как

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

В результате мы имеем кубическую кривую Безье - В-сплайн. Соответствующие базисные функции изображены на рис. 5-27b. На рис. 5-37 приведен еще один пример открытых базисных функций.

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

Векторы могут быть периодическими или открытыми, например

,

На рис. 5.38b-е показаны примеры неравномерных базисных функций В-сплайна порядка . У соответствующих узловых векторов на концах находится по совмещенных одинаковых значений. Для сравнения на рис. 5-38а приведены базисные функции для открытого равномерного вектора. Отметим, что на рис. 5-38а и b функции симметричны, а также что у неравномерных базисов симметрия нарушается: 5-38с-е. Кроме того, при совмещенных узловых значениях у одной из функций появляется излом. На рис. 5-38d и е видно, что положение излома зависит от расположения совмещенного значения в узловом векторе.

Формула Кокса-де Бура (5-84) для расчета базисных функций В-сплайна рекурсивна, поэтому функция порядка зависит от базисных функций более низкого порядка вплоть до 1. Пусть дана базисная функция . Тогда эту зависимость можно выразить в виде треугольника

.

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

.

Рассмотрим пример расчета базисных функций.

Пример 5-9 Расчет периодических базисных функций Найти четыре базисные функции , , третьего порядка . Число в данном случае равно 4. Зависимость базисных функций для изображается следующей диаграммой: Рис. 5-38 Функции неравномерного базиса для , . Обратные зависимости для : . Каков диапазон узлового вектора, необходимый для этого расчета? Из уравнения (5-84) следует, что для вычисления необходимы узловые величины и , а для - и , т. е. необходимы значения от 0 до . Отсюда количество узловых значений равно . Узловой вектор для заданных периодических функций: где . Диапазон параметра . Используя уравнение (5-84) и приведенные выше диаграммы, получим базисные функции для различных значений параметра: ; ; , , ; ; , ; ; , , ; ; ; , ; ; , ; ; ; , ; ; , , ; ; ; , , ; , . Знак в определении приводит к тому, что при все функции равны 0. Результаты показаны на рис. 5-36 и 5-39с. Отметим, что каждая из базисных функций является кусочной параболической (квадратичной) кривой. Параболические сегменты на интервалах , , объединяются и составляют базисные функции . Каждая функция представляет собой параллельный перенос другой.

Пример 5-9 показывает, как построить базис по функциям базиса более низкого порядка. На рис. 5-39а изображены функции первого порядка из примера 5-9, на рис. 5-39b - второго порядка, и на рис. 5-39с - третьего порядка. Обратим внимание на то, как растягивается диапазон ненулевых значений функций с увеличением их порядка. Говорят, что функция базиса обеспечивает поддержку на интервале от до .

Внимательно рассматривая рис. 5-36, можно заметить важное свойство функций равномерного базиса. Из уравнения (5-85) известно, что для любого значения параметра . Отсюда следует, что все множество периодических базисных функций для определено только в диапазоне . За его границами . Для равномерного начинающегося с 0 узлового вектора с целыми приращениями пригодный диапазон параметра: . Для более общих или нормализованных векторов сокращение диапазона параметра соответствует потере интервала узловых значений на каждом конце вектора.

Рис. 5-39 Построение периодических базисных функций . (а) ; (b) ; (c).

Пример 5-10 Расчет открытого равномерного базиса Найти четыре базисные функции , , третьего порядка . Напомним, что формально открытый узловой вектор с целыми интервалами между внутренними узловыми значениями определяется как Диапазон изменения параметра , т. е. от 0 до максимального узлового значения. Как и в примере 5-9, количество узловых значений равно . Если брать узловые значения, то вектор для данного примера примет вид где . Параметр изменяется от 0 до 2. Пользуясь уравнениями (5-84) и диаграммами зависимости, получаем функции базиса для различных диапазонов параметра: ; ; ; ; , , ; ; ; , . Эти результаты приведены на рис. 5-40.

Сравнивая результаты примера 5-10 (рис. 5-40) и 5-9 (рис. 5-39), мы видим, что они существенно различаются для периодического и открытого равномерного узловых векторов. В частности, отметим, что у открытых равномерных узловых векторов на всем диапазоне изменения параметра определен полный набор базисных функций; т.е. для всех . У периодического вектора диапазон параметра уменьшается.

Рис. 5-40 Построение открытых базисных функций . (а) ; (b) ; (с) .

Пример 5-11 Расчет неравномерных базисных функций Найти пять базисных функций , , третьего порядка с узловым вектором , содержащим повторяющееся внутреннее узловое значение. Из уравнений (5-84) и диаграмм зависимости ; ; В частности, как следствие повторения узлового значения, для всех . ; ; , , ; ; , . Результат приведен на рис.5-38d. Заметим, что для всех значений имеем . Например, для , Аналогично для

Рис. 5-41 Зависимость формы В-сплайна от его порядка.

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

Гибкость базиса В-сплайна позволяет воздействовать на форму кривой разными способами:

Изменяя тип узлового вектора и базиса: периодический равномерный, открытый равномерный и неравномерный.

Меняя порядок базисных функций.

Меняя количество и расположение вершин определяющего многоугольника.

Используя повторяющиеся вершины.

Используя повторяющиеся узловые значения в узловых векторах.

Рассмотрим эти способы сначала для открытых В-сплайнов, затем для равномерных периодических и неравномерных В-сплайнов.

Открытый В-сплайн по своим свойствам во многом аналогичен кривым Безье. Как уже отмечалось, если порядок В-сплайна равен количеству вершин определяющего многоугольника, то базис В-сплайна сводится к базису Бернштейна, а сам В-сплайн становится кривой Безье. У открытого В-сплайна любого порядка первая и последняя точки кривой совпадают с соответственными вершинами многоугольника. Далее, наклон кривой в первой и последней вершинах многоугольника равен наклону соответственных сторон многоугольника.

На рис. 5-41 изображены три открытых В-сплайна различного порядка, заданные одним набором из четырех вершин. Кривая четвертого порядка - это кривая Безье - один кубический полиномиальный сегмент. Кривая третьего порядка состоит из двух параболических сегментов, соединяющихся в центре второго отрезка с непрерывностью . Кривая второго порядка совпадает с определяющим многоугольником. Она состоит из трех линейных сегментов, соединяющихся во второй и третьей вершинах с непрерывностью . Угол наклона на концах, заданный наклоном сторон многоугольника, одинаков для всех трех кривых.

Рис. 5-42 Влияние кратности вершины на форму В-сплайна, .

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

На рис. 5-42 изображается влияние повторяющихся или совпадающих вершин. Все В-сплайны имеют порядок . Нижняя кривая определена четырьмя вершинами с узловым вектором . У средней кривой пять определяющих вершин, причем две из них повторяются во второй вершине многоугольника . Узловой вектор - . Верхняя кривая определена шестью вершинами с тремя повторяющимися в точке . Узловой вектор - . Соответствующие многоугольники для трех кривых соответственно таковы: ; и .

Нижняя кривая состоит из единственного кубического сегмента. Средняя кривая состоит из двух сегментов, соединенных между и . Верхняя кривая состоит из трех сегментов: первый от до , второй от до середины между и , третий от этой точки до . Обратим внимание на то, что с увеличением кратности вершины кривая все ближе подходит к . Когда кратность достигается , возникает острый угол в соответствии со свойством выпуклой оболочки В-сплайнов. При внимательном изучении рис. 5-42 можно заметить, что на обеих сторонах кратной вершины имеется линейный участок.

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

Наконец, заметим, что у всех кривых наклон в концах одинаков.

На рис. 5-43 показаны три В-сплайна четвертого порядка. Каждый определяющий многоугольник состоит из восьми вершин. Кривые отличаются тем, что точка передвигается в и . Передвижение точки воздействует на кривую только локально: изменяются лишь сегменты, отвечающие отрезкам , и , .

Рис. 5-43 Локальная коррекция В-сплайна.

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

Для иллюстрации рассмотрим пример.

Пример 5-12 Расчет открытого В-сплайна Рассмотрим многоугольник из примера 5-7: , , , . Найти В-сплайн второго и четвертого порядка. Для открытый узловой вектор , где , , …, . Параметр изменяется от 0 до 3. Кривая состоит из трех линейных сегментов. Для функции базиса имеют вид: ; ; , , ; ; , . Для каждого из этих интервалов В каждом случае в результате получается уравнение прямой для ребра многоугольника, т.е. кривая совпадает с многоугольником. Последняя точка кривой требует особого внимания. Из-за того что интервал в уравнении (5-84а) открыт справа, все базисные функции при равны нулю. Следовательно, последняя точка многоугольника формально не лежит на кривой. Фактически же это не так. Рассмотрим , где - бесконечно малая величина. Если стремится к нулю, то в пределе последние точки кривой и многоугольника совпадают. На практике либо явно добавляют последнюю точку в описание кривой, либо определяют . Для порядок кривой совпадает с количеством вершин определяющего многоугольника, поэтому В-сплайн сводится к кривой Безье. Узловой вектор с . Функции базиса таковы: ; ; , , ; ; ; ; Из уравнения (5-83) получаем параметрический В-сплайн Итак, при Сравнивая с примером 5-7, видим, что результаты одинаковы. Получившаяся кривая изображена на рис. 5-28.

Рис. 5-44 Зависимость формы периодического В-сплайна от его порядка.

Теперь займемся периодическими В-сплайнами. На рис. 5-44 показаны три периодических В-сплайна разного порядка. Все кривые определены теми же вершинами, что и для открытых В-сплайнов на рис. 5-41. Для В-сплайн опять совпадает с определяющим многоугольником. Отметим, однако, что у периодического сплайна при первая и последняя точки на кривой не совпадают с первой и последней точками многоугольника. Наклон в первой и последней точках также может отличаться от наклона соответствующих сторон многоугольника. Для В-сплайн начинается в середине первого ребра и оканчивается в середине последнего, как отмечено стрелками. Это происходит из-за сокращения диапазона параметра для базисных функций периодического В-сплайна. Для периодический узловой вектор - с диапазоном параметра . Для периодический узловой вектор - с диапазоном параметра . Для периодический узловой вектор - с диапазоном параметра .

Сравнение полученных результатов и результатов для открытых узловых векторов на рис. 5-41 показывает, что кривую можно определить на полном диапазоне параметра, задавая кратные узловые значения на концах векторов. При этом кривая растягивается к концам многоугольника.

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

Рис. 5-45 показывает, что воздействие кратных вершин одинаково для периодических и открытых В-сплайнов. Область вокруг точки изображена более подробно.

Рис. 5-45 Влияние кратности вершины на форму периодического В-сплайна, .

Пример 5-13 Расчет периодического В-сплайна Снова рассмотрим многоугольник на рис. 5-44. Вершины многоугольника - , , , . Найти периодический В-сплайн четвертого порядка , заданный этим многоугольником. Для на диапазоне узловой вектор для периодических базисных функций - . Для этого диапазона базисные функции первого порядка (см. уравнение (5-84а), : Из уравнения (5-84b) получаем базисные функции высших порядков: ; ; , , ; ; ; , , ; . , . Точка на В-сплайне при имеет вид Полная кривая приведена на рис. 5-44.

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

Первая вершина совпадает с последней. В-сплайн получается незамкнутым из-за ограниченного диапазона параметра. Здесь периодический равномерный узловой вектор: с диапазоном параметра .

Повторяя всего вершин в начале и/или конце замкнутого многоугольника, получаем замкнутый периодический В-сплайн. (Далее в этом разделе обсуждается другой метод в матричной формулировке.) Результат показан на рис. 5-46b, где определяющий многоугольник - с периодическим равномерным узловым вектором и диапазоном параметра . Альтернативные многоугольники или приводят к тому же самому результату.

На рис. 5-47 одна вершина многоугольника передвинута на другое место. Ее воздействие распространяется на сегменты кривой, соответствующие ребрам многоугольника по обе стороны от сдвинутой точки. На рис. 5-48 показан эффект кратной вершины . Область вокруг изображена более подробно. Опять отметим, что ее влияние распространяется также только на сегмента кривой по обе стороны кратной вершины.

Теперь рассмотрим неравномерные В-сплайны. На рис. 5-49 кривая изменяется под воздействием кратных внутренних узловых значений.

Рис. 5-46 Замкнутый периодический В-сплайн. (а) определяющий многоугольник; (b) - определяющий многоугольник.

Рис. 5-47 Изменение формы В-сплайна после перемещения одной вершины многоугольника.

Верхняя кривая третьего порядка рассчитана для открытого узлового вектора . Базисные функции для этой кривой изображены на рис. 5-38а. Нижняя кривая третьего порядка построена с неравномерным узловым вектором . Ее базисные функции - на рис. 5-38d.

Рис. 5-48 Влияние кратных вершин на форму замкнутого периодического В-сплайна.

Такая же кривая получается с неравномерным узловым вектором и базисными функциями на рис. 5-38е.

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

Открытые неравномерные В-сплайны третьего порядка на рис. 5-50 построены с помощью узлового вектора с внутренними значениями, пропорциональными длинам ребер между вершинами многоугольника. Узловой вектор имеет вид

, ,

, , (5-86)

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

Рис. 5-49 Неравномерные В-сплайны, . (а) ; (b) .

Рис. 5-50 Сравнение открытых неравномерных В-сплайнов: (a) равномерный узловой вектор; (b) неравномерный узловой вектор, пропорциональный длинам хорд; (c) неравномерный узловой вектор, пропорциональный длинам хорд, с двойной вершиной в .

Для сравнения приведена кривая с открытым равномерным вектором, а также кривая с парой совпадающих вершин .

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

Пример 5-14 Неравномерный В-сплайн Найти открытый В-сплайн третьего порядка, заданный точками , , , , с помощью неравномерного вектора с узловыми значениями, пропорциональными длинам ребер многоугольника. Сначала найдем длины хорд Их суммарная длина Из уравнения (5-86) получаем внутренние узловые значения . Узловой вектор имеет вид где , , …, . Диапазон параметра . Кривая состоит из трех параболических сегментов. Функции базиса для таковы ; ; , , ; ; ; , . Функции базиса для :, на единицу меньше количества вершин. Заметим, что уравнение (5-87) содержит только: и равны нулю. Используя это свойство, можно значительно сэкономить вычисления., кривой каждая базисная функция - кусочно кубическая, первая производная - кусочно параболическая, вторая производная - кусочно линейная. Третья производная разрывна и состоит из различных констант., .

Машинное обучение

• Tutorial

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


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

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

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

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

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

где - B-сплайн k+1-го порядка :

Вот как, например, выглядит базис на сетке из g = 9 узлов, равномерно распределенных на интервале :


Сходу разобраться в построении сплайнов через B-сплайны очень сложно. Больше информации можно найти .

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

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

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

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

Пусть также

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

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

Сглаживание

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


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

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

Давайте рассмотрим ошибку
Здесь q - вес функции, влияющей на сглаживание, и

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

где

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


Но при q = 1e-7 в данном примере уже достигается достаточное сглаживание.

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

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


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

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

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

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

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

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

и

Объяснение алгоритма:

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

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

Константы c 0 и c 1 находятся из условий и .

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

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

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

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


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

Теги:

• сплайны
• оптимизация
• python3

Добавить метки

Построение интерполяционных сплайнов в B-форме

Интерполяция сплайном в B -форме позволяет просто управлять гладкостью сплайна. Сплайн в B -форме определяется последовательностью узлов {t i } , среди которых могут быть повторяющиеся. Число повторений узла определяет гладкость сплайна в точках разрыва . Во внутренних точках разрыва должно выполняться следующее правило: число повторений узлов в последовательности + число условий непрерывности в точке разрыва = порядку сплайна. Отметим, что полюса интерполяции (т.е. точки, в которых сплайн совпадает с интерполируемой функцией) вовсе не обязаны совпадать с точками разрыва , хотя произвольно они выбираться не могут!

Интерполирование при помощи функции spapi

Рассмотрим простой пример. Проинтерполируем функцию x sin x на отрезке сплайном f(x) четвертого порядка (k=4) , выбрав семь полюсов интерполяции

.
Обычно первые k узлов совпадают с первым полюсом, а последние k узлов - с последним. Возьмем следующую последовательность из одиннадцати узлов (число узлов должно равняться числу полюсов интерполяции + порядок сплайна, в данном случае 11=7+4)

.
и применим функцию spapi: tau = ; % полюса интерполяции y = tau.*sin(tau); % значения интерполируемой функции в полюсах интерполяции t = ; % последовательность узлов сплайна sp = spapi(t, tau, y); % конструирование сплайна figure fnplt(sp) % построение графика сплайна hold on plot(tau, y, "ro") % отображение данных маркерами Получаем следующий график (см. рис. 1), где сплайн f(x) нарисован линией, а маркеры соответствуют значениям интерполируемой функции x sin x в полюсах интерполяции :

Рис. 1. Интерполяция сплайном f(x) в B -форме при помощи spapi.
Убедимся, что внутренние узлы являются точками разрыва. Их кратность равна единице, следовательно, число условий непрерывности (порядок сплайна минус кратность узла, т.е 4-1) в каждой из них равно трем и, следовательно, должна быть непрерывна сама сплайн-функция, ее первая и вторая производные. Для этого построим графики первой, второй и третьей производных сплайна f(x) , вычислив значения производных сплайна при помощи функции fnder. На графики поместим вертикальные линии в точках разрыва (см. рис. 2): figure subplot(3,1,1) fnplt(fnder(sp,1)) % находим первую производную сплайна и строим ее график title("1st derivative of the spline") hold on lims=axis; ymin=lims(3); ymax=lims(4); % узнаем пределы осей plot(, ,"k") % рисуем линии в точках разрыва subplot(3,1,2) fnplt(fnder(sp,2)) % находим вторую производную сплайна и строим ее график title("2nd derivative of the spline") hold on lims=axis; ymin=lims(3); ymax=lims(4); % узнаем пределы осей plot(, ,"k") % рисуем линии в точках разрыва subplot(3,1,3) fnplt(fnder(sp,3)) % находим третью производную сплайна и строим ее график title("3rd derivative of the spline") % рисуем линии в точках разрыва hold on lims=axis; ymin=lims(3); ymax=lims(4); % узнаем пределы осей plot(, ,"k")
Рис. 2. Графики первой, второй и третьей производной сплайна f(x) в B -форме.
Итак, узлы простой кратности действительно являются точками разрыва третьей производной сплайна.
Если бы требовалось приблизить функцию сплайном, у которого вторая производная в точке 1.7 была бы разрывной, то необходимо было повторить 1.7 дважды в последовательности узлов:

.
t = sp = spapi(t, tau, y); Графики сплайна и его производных строится аналогично, из графиков производных видно (см. рис. 3), что вторая производная сплайна имеет разрыв в узле 1.7:

Рис. 3. Графики первой, второй и третьей производных сплайна в B -форме для другой последовательности узлов (узел 1.7 имеет кратность 2).

B-форма, базисные сплайны

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

Сплайн f(x) должен совпадать с функцией функции x sin x в полюсах интерполяции

Tau = ; % полюса интерполяции y = tau.*sin(tau); % значения интерполируемой функции в полюсах интерполяции t = ; % последовательность узлов сплайна sp = spapi(t, tau, y) % конструирование сплайна Обратимся к содержимому структуры sp, в которую записана информация о B -форме сплайна f(x) : sp = form: "B-" knots: coefs: number: 7 order: 4 dim: 1 Самыми важными для нас сейчас является поля number и coefs - оказывается B -форма состоит из семи сплайнов (sp.number=7), назовем их B 1 (x) , B 2 (x) , ... , B 7 (x) . Сплайн-функция f(x) является суммой этих сплайнов, умноженных на соответствующие коэффициенты C i :

Массив коэффициентов C i записан в sp.coefs. Что же это за сплайны B 1 (x) , B 2 (x) , ... , B 7 (x) ? Оказывается каждый из них построен по своей последовательности из пяти узлов, указанной в таблице:

т.е. сплайн B i (x) построен по узлам t i , t i+1 , ..., t i+k . Определение B -сплайна через его узлы мы рассмотрим позже. Сейчас нам понадобится функция spmak, входящая в состав Spline Toolbox, которая конструирует B -сплайн по заданной последовательности узлов. Для построения, например, сплайна B i (x) следует обратиться к ней так: B1 = spmak(, 1), или задать узлы при помощи сечения массива t: B1 = spmak(, 1). Результатом является структура B1 с информацией о сплайне. Создадим сплайны B 1 (x) , B 2 (x) , ... , B 7 (x) , запишем их в структуры B1, B2, ..., B7, соответственно, и построим их графики разыми цветами (см. рис. 4) при помощи функции fnplt (отдельные структуры для сплайнов взяты только для наглядности изложения, можно было ввести массив структур и все сделать в циклах): % Создание сплайнов B1, B2, ... , B7 B1 = spmak(t(1:5), 1); B2 = spmak(t(2:6), 1); B3 = spmak(t(3:7), 1); B4 = spmak(t(4:8), 1); B5 = spmak(t(5:9), 1); B6 = spmak(t(6:10), 1); B7 = spmak(t(7:11), 1); % Построение графиков сплайнов B1, B2, ... , B7 figure fnplt(B1, "r"); hold on; fnplt(B2, "g"); fnplt(B3, "b"); fnplt(B4, "c"); fnplt(B5, "m"); fnplt(B6, "y"); fnplt(B7, "k"); % Подписываем координаты узлов сплайна set(gca, "XTick", t(4:8), "XTickLabel", num2str(t(4:8)"), "XGrid", "on") % Добавление легенды str = {"spline B_1(x)"; "spline B_2(x)"; "spline B_3(x)";... "spline B_4(x)"; "spline B_5(x)"; "spline B_6(x)"; "spline B_7(x)"}; legend(str, -1)

Рис. 4. Графики B-сплайнов, образующих сплайн f(x) .
Вне интервала (t i , t i+k) B -сплайн B i (x) тождественно равен нулю, в чем несложно убедиться, построив их графики (см. рис. 5), например: figure % строим график B1 на отрезке [-1, 4] subplot(2, 1 ,1) fnplt(B1, [-1 4], "r") title("spline B_1") set(gca,"XTick", , "XTickLabel", num2str("), "XGrid", "on") % строим график B2 на отрезке [-1, 4] subplot(2, 1, 2) fnplt(B2, [-1 4], "g") title("spline B_2") set(gca, "XTick", , "XTickLabel", num2str("), "XGrid", "on")

Рис. 5. Графики сплайнов B 1 (x) и B 2 (x) на отрезке [-1, 4].
Обратите внимание, что сплайн B 1 (x) имеет разрыв в точке 0 потому, что его среди его узлов t 1 =0 , t 2 =0 , t 3 =0 , t 4 =0 , t 5 =1.2 узел со значением 0 встречается 4 раза. Так как сам сплайн 4-го порядка, то число условий непрерывности равно нулю (4-4=0), т.е. нет непрерывности не только производных, но и самого сплайна.
Сумма сплайнов B 1 (x) , B 2 (x) , ... , B 7 (x) , умноженных на соответствующие коэффициенты (которые возвращаются в sp.coefs при вызове sp = spapi(t, tau, y)),

является сплайном в B-форме, ровно тем, который возвращает функция spapi. Это несложно проверить, вычислив данную сумму сплайнов. Для нахождения значений каждого из сплайнов B 1 (x) , B 2 (x) , ... , B 7 (x) в точках отрезка с шагом 0.05 применим функцию fnval, значения суммы в этих точках запишем в вектор f и построим график зависимости суммы от x на отрезке : % записываем коэффициенты сплайна в B-форме в вектор c c = sp.coefs; % вычисляем значения сплайна на отрезке с шагом 0.05 x = 0:0.05:3; f = c(1)*fnval(B1, x) + c(2)*fnval(B2, x) + c(3)*fnval(B3, x) + c(4)*fnval(B4, x) + ... c(5)*fnval(B5, x)+c(6)*fnval(B6, x) + c(7)*fnval(B7, x); % строим графики сплайна и исходных данных figure plot(x, f) hold on plot(tau, y, "ro") Построенный график совпадает с графиком сплайна, полученного ранее при непосредственном использовании функции spapi для интерполяции функции x sin x , поэтому график сплайна здесь не приводится (см. рис. 1).

Выбор узлов B-сплайнов

Как мы уже замечали, последовательность узлов B -формы сплайна не может выбираться независимо от последовательности полюсов интерполяции. Во-первых, число узлов должно равняться числу полюсов + порядок сплайна. Во-вторых, для существования и единственности сплайна требуется выполнение еще условия теоремы Шенберга-Уитни.
Начнем с примера. Попробуем взять такую последовательность узлов:

А полюса интерполяции оставим теми же самыми

При выполнении команд tau = [ 0 0.5 1.1 1.8 2.2 2.8 3]; y = tau.*sin(tau); t = sp = spapi(t, tau, y); получим ошибку и сообщение о том, что матрица является вырожденной. Дело в том, что при интерполяции B -сплайнами для нахождения коэффициентов B -формы необходимо решить некоторую систему линейных алгебраических уравнений, матрица которой в данном случае оказывается вырожденной. Для однозначной разрешимости этой системы линейных уравнений, а следовательно, и для однозначного определения интерполяционного сплайна все элементы последовательности полюсов интерполяции и узлов должны удовлетворять неравенствам:

В которых строгое равенство допускается лишь для крайних узлов (разумеется, имеется ввиду, что последовательность узлов упорядочена по возрастанию).
В данном случае уже не выполняется условие , так как Конечно, для построения подходящей последовательности узлов сплайна не требуется подбирать их вручную. В Spline Toolbox входит функция aptknt , которая по заданной последовательности полюсов интерполяции tau возвращает последовательность t узлов, обеспечивающую существование и единственность сплайна заданного порядка k. Обращение к функции aptknt выглядит следующим образом: t = aptknt(tau, k) Для интерполяции B -сплайнами 4-го порядка функции в нашем примере подошла бы такая последовательность команд: tau = [ 0 0.5 1.1 1.8 2.2 2.8 3]; y = tau.*sin(tau); t = aptknt(tau, 4) sp = spapi(t, tau, y); figure fnplt(sp) % построение графика сплайна hold on plot(tau, y, "ro") % отображение данных маркерами Более того, интерфейс функции spapi позволяет вообще не заботится об узлах сплайна и указать только степень, полюса интерполяции и значения функции в них: sp = spapi(4,tau,y); Результат будет тот же самый, поскольку при таком способе обращения к функции spapi она вызывает функцию aptknt для конструирования допустимой последовательности узлов сплайна, удовлетворяющей условиям теоремы Шенберга-Уитни. Функция aptknt возвращает такую последовательность узлов:

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

для i=1,2,...,n-k . Предполагается, что последовательность полюсов интерполяции содержит не менее элементов, является неубывающей и для всех i . Эти узлы удовлетворяют условиям теоремы Шенберга-Уитни и обеспечивают существование и единственность интерполяционного сплайна.
Вообще говоря, выбор хороших полюсов и узлов интерполяции сплайна (оптимальных в некотором смысле) - отдельная задача. В Spline Toolbox имеется специальная функция optknt , которая так же, как и aptknt, по заданным полюсам интерполяции и заданному порядку сплайна возвращает последовательность его узлов, удовлетворяющую условию теоремы Шенберга-Уитни. Узлы сплайна, которые строит функция optknt, приводят к наилучшему приближению сплайном.
Обратная задача - выбор в некотором смысле оптимальных полюсов для интерполяции сплайном заданного порядка с заданными узлами - решается при помощи функции chbpnt . Хорошее приближение также дают полюса интерполяции, получаемые путем усреднения узлов сплайна. Для сплайна порядка k с последовательностью узлов такие полюса определяются по формуле:

Для их вычисления в Spline Toolbox предназначена функция aveknt. Обсуждению проблемы выбора полюсов интерполяции и узлов сплайна будет посвящен специальный раздел.

Кроме перечисленных функций, в состав Spline Toolbox входят несколько сервисных функций для работы с узлами и полюсами интерполяции.

Функция augknt позволяет просто создать последовательность узлов сплайна по заданным точкам разрыва так, что в каждой точке разрыва узел будет иметь заданную кратность, например >> t = augknt([-2 -1 0 1 2], 4) возвращает t = -2 -2 -2 -2 -1 0 1 2 2 2 2 Если требуется задать кратность узла в каждой внутренней точке отдельно, то следует указать вектор кратностей в качестве третьего входного аргумента: >> t = augknt([-2 -1 1 2], 4, ) t = -2 -2 -2 -2 -1 -1 1 1 1 2 2 2 2 Для конструирования последовательности узлов с заданными кратностями подходит также функция brk2knt, которая требует задания двух входных аргументов - массива точек разрыва и массива их кратностей такой же длины: >> knots = brk2knt([-1 0 1], ) knots = -1 -1 -1 0 0 1 1 1 1 Обратная задача - получение точек разрыва и соответствующих кратностей узлов решается при помощи функции knt2brk: >>knots = [-1 -1 -1 0 0 1 1 1 1] >> = knt2brk(knots) breaks = -1 0 1 mults = 3 2 4

Интерполирование B-сплайнами с заданными значениями производных

Функция spapi позволяет конструировать такие B -сплайны, которые принимают заданные значения вместе со своими производными в полюсах интерполяции. Для этого при вызове функции spapi: sp = spapi(t, tau, y), где t - массив узлов, или sp = spapi(k, tau, y), где k - порядок сплайна следует повторить координату полюса в массиве tau и, соответственно, в векторе y задать значения не только сплайна, но и его производной в этом полюсе интерполяции. Для задания значения сплайна и его p -ой производной в некотором полюсе , он должен встречаться p+1 раз в массиве tau, тогда соответствующие значения массива y будут восприниматься как значения сплайна, его первой, второй и т. д., ..., p -ой производной.
Например, для интерполяции кубическим сплайном некоторой функции, заданной таблицей:

подойдет такая последовательность команд (см. результат на рис. 6, на котором приведен график сплайна, его первой производной и маркерами отмечены заданные значения в полюсах интерполяции): x = ; y = ; dy = sp = spapi(4, , ) figure subplot(2, 1, 1) fnplt(sp) title("cubic spline") hold on plot(x, y, "ro") subplot(2,1, 2) fnplt(fnder(sp)) title("its first derivative") hold on plot(x, dy, "ro")
Рис. 6. Интерполяция B-сплайном по заданным значениям функции и ее производной.
Значения производных интерполируемой функции могут быть заданы не во всех полюсах интерполяции. Рассмотрим пример, в котором требуется интерполировать функцию y=e -x sin3x на отрезке в шести полюсах: так, чтобы значения первой и второй производных сплайна f(x) совпадали со значениями y и в трех полюсах . Для этого задаем массив полюсов интерполяции tau, inline-функцию для вычисления и записываем значения ее первой и второй производных в полюсах в массивы y1 и y2, соответственно: fun = inline("exp(-x).*sin(3*x)"); tau = 0:5; y = fun(tau); y1 = exp(-tau(1:2:5)).*(3*cos(3*tau(1:2:5))-sin(3*tau(1:2:5))); y2 = -exp(-tau(1:2:5)).*(8*sin(3*tau(1:2:5))+6*cos(3*tau(1:2:5))); Далее обращаемся к функции spapi: f = spapi(4, , ); Сравним исходную функцию, полученный сплайн и еще два сплайна: g(x) (который был построен с учетом значений производных ) и сплайн h(x) (который просто совпадает с функцией f(x) в полюсах интерполяции). Результат представлен на рис. 7. g = spapi(4, , ); h = spapi(4, tau, y); figure fplot(fun, , "b") hold on fnplt(f, "m") fnplt(g, "g") fnplt(h, "k") legend("y(x)", "f(x)", "g(x)", "h(x)")
Рис. 7. Интерполяция функции y=e -x sin3x с учетом и без учета значений производных.