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

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

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

- (от англ. spline, от spline гибкое лекало, полоса металла, используемая для черчения кривых линий) функция, область определения которой разбита на конечное число отрезков, на каждом из которых сплайн совпадает с некоторым… … Википедия

Сущ., кол во синонимов: 1 функция (49) Словарь синонимов ASIS. В.Н. Тришин. 2013 … Словарь синонимов

Сплайн-функция - кусочно гладкая функция, используемая для выравнивания временных рядов. Применение С. ф. вместо обычных функций тренда эффективно, когда внутри анализируемого периода меняется тенденция, направление ряда. С. ф. помогает… … Экономико-математический словарь

сплайн - (Spline) Математическая кривая, плавно соединяющая отдельные точки. Применяется для изображения контуров знаков [граница изображения знака]. См. также кривые Безье [метод описания веторных кривых] … Шрифтовая терминология

Кубический эрмитов сплайн сплайн, построенный из кубических полиномов с использованием эрмитовой интерполяции, в соответствии с которой интерполируемая функция задается не только своими значениями в n точках, но и ее первыми производными. Для… … Википедия

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

Интерполирование посредством сплайнов, т. е. построение интерполяционного сплайна (и. с.), принимающего в заданных точках {xi}заданные значения {f(xi)}, i=0, 1, . . ., n. Обычно и. с. удовлетворяют дополнительным условиям в концевых точках. Так,… … Математическая энциклопедия

Функция определенная на отрезке ,совпадающая на частичных отрезках [ х i, xi+1], образованных сеткой а=x0Математическая энциклопедия

сплайн - а, ч. Одна з елементарних функцій, включена у сучасний числовий аналіз … Український тлумачний словник

сплайн-апроксимація - іменник жіночого роду …

сплайн-інтерполяція - іменник жіночого роду … Орфографічний словник української мови

Книги

Сплайн-всплески и их реализации , Бурова И.Г.. В данной книге основное внимание уделяется сплайн-всплесковым разложениям первого и второго порядка. Рассмотрены приемы разложения потоков в эрмитовом случае с использованием потока значений…

Любой более или менее сложный чертеж состоит не только из отрезков прямых линий, окружностей и их дуг, но также и из набора кривых линий. Гладкие кривые удобно строить при помощи метода сглаживания кривой типа В-сплайна. B-сплайн это гладкая кривая или, точнее, кривая с непрерывными старшими производными до n-ой, где n порядок сплайна. Заметим, что линия, составленная из В-сплайнов, не будет проходить точно через заданные точки. Подобную кривую составляют из дуг полиномов третьей степени, так как такой полином обеспечивает необходимую непрерывность. Построение линии происходит с помощью итерационной процедуры.

Рассмотрим построение кубического сплайна. Пусть нам даны две соседние точки, через которые проведем кубический полином, но у полинома 4 коэффициента, следовательно нужно еще два дополнительных условия или точки. Для этого прихватим еще две соседние точки. Чем более плавной мы хотим видеть линию, тем сложнее пройти точно через точки. Если в формуле x = q 3 , то достаточно плавности 3.

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

Рассмотрим . Пусть t параметр, по которому пробегаем от точки P i к точке P i+1 . При t = 0 мы находимся в точке P i , при t = 1 в точке P i+1 . Если 0 < t < 1, то мы находимся между P i и P i+1 .

Эта линия в каждой точке имеет систему:

x(t) = ((a 3 t + a 2)t + a 1)t + a 0 , для 0 <= t <= 1

y(t) = ((b 3 t + b 2)t + b 1)t + b 0 , для 0 <= t <= 1 a 3 = (-x i-1 + 3x i - 3x i+1 + x i+2)/6

A 2 = (x i-1 - 2x i + x i+1)/2

A 1 = (-x i-1 + x i+1)/2

A 0 = (x i-1 + 4x i + x i+1)/6

Точки b 3 - b 0 расписывают так же, но вместо x подставляют у. Между P i и P i+1 точки а и b не меняются. Если после последней точки указать первую точку, то система замкнет контур.

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

Следствия

Сглаживание B-сплайнами

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

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

Можно строить достаточно гладкие кривые и поверхности с использованием полиномов. Допустим, что мы хотим построить поверхность в виде графика функции z = z(x, y). Линия y = const на этой поверхности будет представлена линией z = z(x), она будет проходить через последовательность точек (x 0 , z 0), ..., (x i , z i), ..., (x n , z n) с x 0 < ... < x i < ... < x n . Наша цель провести через эти точки составную кривую f(x), имеющую следующие свойства:

на каждом подынтервале x i-1 <= x <= x i , i = 1, 2, ..., n функция f(x) является кубическим полиномом;
ее первые и вторые производные непрерывны в узлах.

Полученная гладкая кривая называется кубическим сплайном. Термин «сплайн» возник по аналогии: это название чертежного инструмента тонкой металлической линейки, которая может изгибаться так, чтобы проходить через заданные точки. Физически такая кривая минимизирует энергию внутренних напряжений. Математически имеет минимальную среднеквадратичную кривизну, то есть она наиболее гладкая. Сплайны имеют много приложений в конструировании криволинейных форм. Однако они имеют и некоторые ограничения:

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

Первое ограничение можно устранить с помощью B-сплайна. Общая форма полученной в этом случае кривой показана на .

На этом рисунке сплайн продолжен от его конечных точек x i-4 , x i прямыми линиями, идущими вдоль оси x. В результате получается кубический сплайн на любом числе отрезков, но он не равен нулю только на четырех из них. Такая функция называется B-сплайном (или фундаментальным сплайном) четвертого порядка (или третьей степени). Про него говорят, что он имеет минимальный носитель (носитель это число отрезков, на которых сплайн отличен от нуля).

Заметим, что кубический B-сплайн полностью определяется множеством узлов, на которых он определен, и только одной заданной величиной z. В более общем виде B-сплайн M mi (x) порядка m (или степени m - 1) на данном множестве узлов везде равен нулю, кроме m последовательных отрезков x i-m < x < x i . Опять-таки M mi (x) определяется множеством узлов и одной величиной z. Принято исключать последнюю степень свободы и фиксировать амплитуду B-сплайна некоторым стандартным образом.

Часто удобно для вычислений использовать нормализованный B-сплайн N mi (x), связанный с M mi (x) соотношением N mi (x) = (x i - x i-m)M mi (x).

Любой сплайн порядка m на множестве узлов x 0 , x 1 , ..., x n может быть выражен в виде линейной комбинации B-сплайнов, определенных на том же множестве узлов, расширенном (m - 1) дополнительными узлами на каждом из концов интервала, которые можно выбрать произвольно: x -m+1 , x -m+2 , ..., x -1 и x n+1 , ..., x n+m-1 . Можно построить m + n - 1 последовательных B-сплайнов на расширенном множестве узлов, каждый из которых отличен от нуля на m последовательных отрезках. Поэтому можно записать:
j (x) = S c i * M mi (x),
где j (x) любой сплайн степени (m - 1) на первоначальном множестве узлов и M mi (x) есть B-сплайн на расширенном множестве узлов, отличный от нуля при x i-m < x < x i ; c i суть числовые коэффициенты; суммирование ведется по i = 1, ..., m + n - 1.

Если имеется множество векторов r 0 , r 1 , ..., r n , то можно использовать их: r(u) = S r i * N 4, i+1 (u) (суммирование ведется по i = 0, ..., n). Так как имеется (n + 1) векторных коэффициентов, то необходим набор из (n + 1) B-сплайнов. Последняя формула для 0 <= u <= n - 2 является уравнением кривой, образованной кубическими B-сплайнами.

Свойства

Некоторые простейшие свойства следуют из тождества S N 4, i+1 = 1, 0 <= u <= n - 2, i = 0..n. При u = 0 следует: r(0) = N 42 (0)(r 1 - r 0). Из этого следует, что если r 0 , r 1 , .., r n вершины некоторой замкнутой ломанной, то кривая, построенная на основе B-сплайна, начинается в r 0 и ее касательная в этой точке имеет направление (r 1 - r 0). Аналогичное утверждение верно и для другого конца. Главное преимущество этого сплайна заключается в том, что изменение одной из вершин влечет за собой изменение только четырех отрезков кривой. Далее, мы также можем построить кривую, аппроксимирующую ломанную с любым желаемым числом сторон. Отрезок сплайна всегда лежит в выпуклой оболочке:

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

Имеется еще 2 полезных факта:

кривая проходит вблизи средней точки каждой стороны, за исключением 1-ой и последней точками;
при k = 2, ..., n - 2 кривая проходит через точки: 1/6r k-1 + 2/3r k + 1/6r k+1 = 2/3r k + 1/3(1/2(r k-1 + r k+1))

Эти точки, как показано на , лежат на 1/3 расстояния от r k на прямой, соединяющей r k с серединой отрезка между r k-1 и r k+1 .

Построение интерполяционных сплайнов в 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 с учетом и без учета значений производных.