![](https://habrastorage.org/files/f0d/403/8c0/f0d4038c02204f1cbb7cba71c5dba9fc.png)
![](https://habrastorage.org/files/3fc/d83/f89/3fcd83f893424620a0d1a78d9b8e1c4d.png)
![](https://habrastorage.org/files/767/93c/7de/76793c7dedd64b16ae3e957bc852163c.png)
Под катом находится алгоритм, раскрывающий, каким образом сплайны позволяют строить подобную красивую регрессию:
![](https://habrastorage.org/files/e3c/e94/ff7/e3ce94ff72bd4305837c9d530b8cb4a4.png)
Основные определения
Функция s(x) на интервале [a, b] называется сплайном степени k на сетке с горизонтальными узлами
![](https://habrastorage.org/files/eab/a4f/26a/eaba4f26aea143d6aed5aaf7944a7baf.png)
- На интервалах
функция s(x) является полиномом k-й степени.
- n-ая производная функции s(x) непрерывна в любой точке [a, b] для любого n = 1,…, k-1.
Заметим, что для построения сплайна нужно для начала задать сетку из горизонтальных узлов. Расположим их таким образом, чтобы внутри интервала (a, b) стояло g узлов, а по краям — k+1:
![](https://habrastorage.org/files/c59/cf0/cda/c59cf0cdaadc4f55965f1cc2cc0ec020.png)
![](https://habrastorage.org/files/9e8/5eb/1ea/9e85eb1ea5264a2e90f9f3ba3a842145.png)
Каждый сплайн в точке
![](https://habrastorage.org/files/1fc/a36/8c6/1fca368c69b04ca38cfd6920d12c7218.png)
![](https://habrastorage.org/files/bf0/3a9/022/bf03a902279e4f5fa11a848b6066487b.png)
где
![](https://habrastorage.org/files/199/957/081/199957081ea943ff90133bc9e3221c35.png)
![](https://habrastorage.org/files/872/584/874/872584874fd641e0b646b529d084dc0c.png)
![](https://habrastorage.org/files/d5a/b7f/882/d5ab7f8826684361aaabd0c691814276.png)
![](https://habrastorage.org/files/e21/e92/8c3/e21e928c3cef4d8b8c2bfb9d74b7184e.png)
Вот как, например, выглядит базис на сетке из g = 9 узлов, равномерно распределенных на интервале [0, 1]:
![](https://habrastorage.org/files/9cb/fc2/6b3/9cbfc26b33a4471cad716a718a464d39.png)
Сходу разобраться в построении сплайнов через B-сплайны очень сложно. Больше информации можно найти здесь.
Аппроксимация с заданными горизонтальными узлами
Итак, мы выяснили что сплайн определяется однозначно узлами и коэффициентами. Допустим, что узлы
![](https://habrastorage.org/files/eab/a4f/26a/eaba4f26aea143d6aed5aaf7944a7baf.png)
![](https://habrastorage.org/files/f0d/403/8c0/f0d4038c02204f1cbb7cba71c5dba9fc.png)
![](https://habrastorage.org/files/3fc/d83/f89/3fcd83f893424620a0d1a78d9b8e1c4d.png)
![](https://habrastorage.org/files/61a/9d2/9df/61a9d29dffdd46d3a6a57bc978d9f1ce.png)
![](https://habrastorage.org/files/a7e/2f5/ce2/a7e2f5ce2b5447988640e14f84690a5f.png)
Для удобства запишем в матричном виде:
![](https://habrastorage.org/files/214/cd7/e5a/214cd7e5a62d44fc852e89e655e4f057.png)
где
![](https://habrastorage.org/files/238/8d2/8b3/2388d28b33424d4bb01f40c62329dca8.png)
![](https://habrastorage.org/files/f5b/f44/23c/f5bf4423cd7a44c787b2e177f5955ca3.png)
Заметим, что матрица E блочно-диагональная. Минимум достигается когда градиент ошибки по коэффициентам будет равен нулю:
![](https://habrastorage.org/files/def/3f9/60f/def3f960f95c404890b97d5bb097ee3c.png)
Зададим оператор
![](https://habrastorage.org/files/c91/ee4/982/c91ee4982a2046beb5aea4b4252b8fb6.png)
![](https://habrastorage.org/files/72e/e09/bcb/72ee09bcb9434b8d995f457d8c960526.png)
![](https://habrastorage.org/files/79b/39f/9e7/79b39f9e774f4b1e9c5388b66c3d31ba.png)
Пусть также
![](https://habrastorage.org/files/e47/5fc/12f/e475fc12fb73485ea25db1167dae04a3.png)
![](https://habrastorage.org/files/c0e/44b/227/c0e44b2279344ecabbaaa3779a651f8f.png)
Тогда вся задача и все предыдущие формулы сводятся к решению простой системы линейных уравнений:
![](https://habrastorage.org/files/6c9/aa4/5dd/6c9aa45ddc144cdc95131c613bf9b46d.png)
где матрица А (2k+1)-диагональная, так как
![](https://habrastorage.org/files/b25/e32/f26/b25e32f26ff7419aba245d38d2f7d2a7.png)
И вот, решая систему, получаем желаемый результат:
![](https://habrastorage.org/files/e3c/e94/ff7/e3ce94ff72bd4305837c9d530b8cb4a4.png)
Сглаживание
Однако, далеко не всегда все так хорошо. При малом количестве данных по отношению к количеству узлов и степени сплайна может возникнуть проблема т.н. сверхподгонки (overfitting). Вот пример «плохого» кубического сплайна, при этом идеально проходящего сквозь данные:
![](https://habrastorage.org/files/940/b7a/eba/940b7aebae8d4602bdb1d3bc670f1870.png)
Окей, кривая уже не такая уж и красивая. Попытаемся уменьшить так называемые колебания сплайна. Для этого мы попробуем «сгладить» его k-ю производную. Другими словами, мы минимизируем разницу между производной слева и производной справа от каждого узла:
![](https://habrastorage.org/files/eb4/dcc/a4b/eb4dcca4b0764a86b60d1610758eaf67.png)
Разложив сплайн в базисную форму, мы получаем:
![](https://habrastorage.org/files/2d5/89e/a39/2d589ea39a1048cf8bb9b908f3a1eab4.png)
![](https://habrastorage.org/files/edd/826/15a/edd82615ac1343978f6149c7e72848bd.png)
Давайте рассмотрим ошибку
![](https://habrastorage.org/files/943/f3a/8bc/943f3a8bc46346dc9d08e7f37202b0c5.png)
Здесь q — вес функции, влияющей на сглаживание, и
![](https://habrastorage.org/files/c63/b60/da6/c63b60da69224b00bec24aeb66d6b2c9.png)
![](https://habrastorage.org/files/ea2/7c4/87b/ea27c487b2f943c1b647054e9e8acce6.png)
Новая система уравнений:
![](https://habrastorage.org/files/b48/686/35a/b4868635a1af4aedb962633290bf6331.png)
где
![](https://habrastorage.org/files/299/b53/260/299b532601c04f3b837210d86e699717.png)
Ранг матрицы B равен g. Она симметричная и, так как q > 0, A + qB будет положительно определенной. Поэтому разложение Холецкого по-прежнему применимо к новой системе уравнений. Однако, матрица B вырожденная и при слишком больших значениях q могут возникнуть численные ошибки.
При совсем маленьком значении q = 1e-9 вид кривой изменяется очень слабо.
![](https://habrastorage.org/files/75b/70f/6fc/75b70f6fc8094d9aa91d46d3677dea74.png)
Но при q = 1e-7 в данном примере уже достигается достаточное сглаживание.
![](https://habrastorage.org/files/cbd/af8/918/cbdaf8918b274f03b854757993dbbe7a.png)
Аппроксимация с неизвестными горизонтальными узлами
Представим теперь, что задача такая же, как и прежде, за исключением того, что мы не знаем как узлы расположены на сетке. На вход кроме данных подается только количество узлов g, интервал [a, b] и степень сплайна k. Попробуем наивно предположить, что лучше всего расположить узлы равномерно на интервале:
![](https://habrastorage.org/files/476/81f/56f/47681f56f4ad48d08f1a168e996bd6d9.png)
Упс. Видимо, необходимо все-таки расположить узлы как-то иначе. Формально, расположим узлы таким образом, чтобы значение ошибки
![](https://habrastorage.org/files/c0a/7c1/6c4/c0a7c16c44bc4f5b99d7530c81882de7.png)
было минимально. Последнее слагаемое играет роль штрафной функции, чтобы узлы не сильно приближались друг к другу:
![](https://habrastorage.org/files/c2f/6f0/fcb/c2f6f0fcbc154a9aba09f386540cb811.png)
Положительный параметр p — вес штрафной функции. Чем больше его значение, тем быстрее узлы будут удаляться друг от друга и стремиться к равномерному расположению.
Для решения данной задачи мы используем метод сопряженных градиентов. Его прелесть заключается в том, что для квадратичной функции он сходится за фиксированное (в данном случае g) количество шагов.
- Инициализируем направление
.
Как рассчитать производную ошибки по узлам?Производная суммы квадратов по узлу:
Для того, чтобы рассчитать влияние положения узла на значения сплайна, нужно рассмотреть B-сплайнына новых узлах
и с новыми коэффициентами
Производная штрафной функции:
На производную функции сглаживания без слез не взглянешь:
- Для j = 0,…, g-1
- Зададим функцию
возвращающую ошибку в зависимости от выбора шага вдоль заданного направления. На этом шаге мы находим оптимальное значение ?*, доставляющее минимум этой функции. Для этого мы решаем задачу одномерной оптимизации. О том, каким образом, будет сказано позже.
- Обновляем значения узлов:
- Обновляем вектор направления:
- Зададим функцию
- Если
и
где ?1 и ?2 — заранее заданные величины, отвечающие за точность работы алгоритма, то выходим. Иначе, обнуляем счетчик и возвращаемся на первый шаг.
Решение задачи одномерной минимизации
Для того, чтобы найти значение
![](https://habrastorage.org/files/950/af0/84b/950af084b8b1434fbb7013548c67eef0.png)
![](https://habrastorage.org/files/fb9/f08/d27/fb9f08d2777b4c258edfb8027db7a86b.png)
мы используем алгоритм, позволяющий сократить количество обращений к оракулу, а именно количество операций аппроксимации с заданными узлами и подсчета функции ошибки. Мы будем использовать нотацию
![](https://habrastorage.org/files/14f/355/58f/14f35558faa54eb8a49465e73a1ad990.png)
- Пусть первая и последняя компоненты вектора направления равны нулю:
. Зададим также максимально возможный шаг вдоль этого направления:
Такой выбор обусловлен тем, что узлы не должны пересекаться.
- Инициируем k = 0 и начальные шаги:
,
,
.
- До тех пор, пока
:
- Задаём
и уменьшаем шаг
- k = k + 1
- Задаём
- Если k > 0, то возвращаем ?* = ?1. Иначе:
- До тех пор, пока
: ?0 = ?1, ?1 = ?2 и
- Возвращаем
, где
— корень уравнения I'(?) = 0 и I(?) — аппроксимация функции ошибки:
где
- До тех пор, пока
Коэффициенты ai и bi могут быть найдены из уравнений
![](https://habrastorage.org/files/1df/b31/ede/1dfb31ede82f4be48811dce41aadb8c2.png)
и
![](https://habrastorage.org/files/6ff/c40/41b/6ffc4041bf6a4ae09a3073f0b2ce6d55.png)
Объяснение алгоритма:
Идея заключается в том, чтобы расставить три точки ?0 < ?1 < ?2 таким образом, чтобы по значениям ошибок, достигаемых в этих точках, можно было построить простую аппроксимирующую функцию и вернуть её минимум. Притом значение ошибки в ?1 должно быть меньше, чем значение ошибки в ?0 и ?2.
Находим начальное приближение ?1 из условия S'(?1)=0, где S(?) — функция вида
![](https://habrastorage.org/files/609/fc2/025/609fc20256e04407a04fa5f08d07b554.png)
Константы c0 и c1 находятся из условий
![](https://habrastorage.org/files/f21/827/296/f21827296c954412a8f48306f71b623d.png)
![](https://habrastorage.org/files/8cf/a02/95c/8cfa0295c09e40159a26a274bdc3e548.png)
Если мы просчитались с начальным приближением, то мы уменьшаем шаг ?1 до тех пор, пока он доставляет большее значение ошибки, чем ?0. Выбор
![](https://habrastorage.org/files/fb0/3e4/090/fb03e409006041698caf2b4973a3befb.png)
![](https://habrastorage.org/files/f1d/a00/0b0/f1da000b0aec4f8cb59975abedad7992.png)
![](https://habrastorage.org/files/5d6/a6d/98b/5d6a6d98b48d41e98ff5ea7bd8d90084.png)
![](https://habrastorage.org/files/c35/cb3/1c9/c35cb31c9c9c4e2ba38f0bb7b1f5d3b7.png)
![](https://habrastorage.org/files/e2d/a52/986/e2da52986f4e4dd6b4ffc18798d92b83.png)
![](https://habrastorage.org/files/5ff/ed2/e5f/5ffed2e5fd8847719863502d52f173b0.png)
Если k > 0, то мы нашли значение ?1, такое что при его выборе значение ошибки будет меньше, чем при выборе ?0 и ?2, и мы возвращаем его в качестве грубого приближения ?*.
Если же наше первоначальное приближение было верным, то мы пытаемся найти шаг ?2, такой что
![](https://habrastorage.org/files/820/1b2/93c/8201b293ce7242e281ca84ff030df83b.png)
Когда найдены все три значения ?0, ?1 и ?2, мы представляем функцию ошибки в виде суммы двух функций, приближающих разность квадратов и функцию штрафа. Функция Q(?) — парабола, чьи коэффициенты могут быть найдены, так как мы знаем её значения в трех точках. Функция R(?) уходит на бесконечность при ?, стремящемся к ?max. Коэффициенты bi также могут быть найдены из системы из трех уравнений. В результате, мы приходим к уравнению, которое может быть приведено к квадратному и легко решено:
![](https://habrastorage.org/files/f94/e98/91d/f94e9891d64542b498220427981fa3fa.png)
И вот, для сравнения, результат оптимально построенного сплайна:
![](https://habrastorage.org/files/d80/2d4/ec3/d802d4ec3af342ffb4aebb2473e29b48.png)
Ну и для тех, кому может пригодиться: реализация на Python.
Комментарии (15)
JIuBeHb
07.12.2016 23:23Приблизительно год назад я решал задачу аппроксимации массива точек в трехмере при помощи сплайн-плоскости. К сожалению, массивы у меня были «плохие» и постоянно получались «колодцы» или «горы» там, где их быть не должно, как в примере автора.
В сети очень много сплайн-алгоритмов и их реализаций на различных языках именно для двумерных случаев (первая половина статьи, например, повторяет методичку, которую нам раздавали на втором курсе), а вот для трёхмерных и, тем более, n-мерных случаев толковых алгоритмов, а тем более реализаций, я лично найти не смог.
Может быть автор или кто-нибудь ещё сможет описать реализацию подобных «гладких» сплайнов для трехмера?The_Freeman
08.12.2016 17:01Посмотрите книгу Curve and Surface Fitting with Splines (Numerical Mathematics and Scientific Computation), автор Paul Dierckx. Основные идеи в статье взяты оттуда.
kanikeev
08.12.2016 01:33Человеческим зрением вижу две кривых на первом рисунке. Вы ничего в критериях не потеряли?
3aicheg
08.12.2016 04:25-7Очередной сэнсэй дорвался до редактора формул и наваял «туториал», ненужный тем, кто умеет и непонятный тем, кто не умеет :)
Refridgerator
08.12.2016 06:26В своё время пытался подобным образом решать аналогичную задачу, но не осилил, заткнулся на выборе узлов. Решил её другим способом — последовательной аппроксимацией функцией a*sin(k*x)+b*cos(k*x)
0serg
08.12.2016 08:59+2Неплохая статья, но предложенный метод решения проблемы наименьших квадратов не самый устойчивый — вычисление AT A возводит число обусловленности проблемы в квадрат. Часто надежнее (но медленнее) для подобных проблем использовать SVD
aso
08.12.2016 09:57На всякий случай скажу очевидную вещь — выбор критериев остаётся делом произвола и не имеет однозначных формальных критериев.
Ну т.е. приведённые «нехорошие примеры» — в реальности, в некоторых случаях — могут оказаться вполне хорошими.
Понимать это может только человек, подбирающий сплайн.
Caduceus
08.12.2016 10:21Я нечто подобное в свое время генетическим алгоритмом делал. Довольно забавно в динамике наблюдать, как кривая подстраивается под точки.
Refridgerator
09.12.2016 07:26+3Чтение математической литературы может быть увлекательнее любого детектива.
Сначала — завязка. Поставлена понятная задача и показано решение. Интересно, открываем.
Функция s(x) на интервале [a, b] называется сплайном степени k на сетке с горизонтальными узлами
Что за g, почему 2k? Интрига. Но уже в следующем предложении интрига раскрывается, заодно показывая несколько иллюстраций, из которых понятно, что «ну это всё просто!». Расслабляемся и читаем дальше.
Строго говоря, они должны доставлять минимум функции
Хмм… Похоже на метод наименьших квадратов. Приятно чувствовать, что какие-то знания у меня есть. И тут вдруг — неожиданный поворот сюжета:
Для удобства запишем в матричном виде:
Что за палочки, что за цифры, нас такому не учили! Страшно. И хотя дальше матрица в привычном представлении присутствует, чувство тревоги не уходит. И вот:
Тогда вся задача и все предыдущие формулы сводятся к решению простой системы линейных уравнений:
Звучит как очевидный факт, но мне же это совсем не очевидно! Расстроившись, я пошёл пить чай.
(спустя некоторое время)
Хм… Минимизировать разницу между производными для обеспечения гладкости — это сильно!
(спустя ещё некоторое время)
Что-то формул всё больше и больше. Посмотрю-ка я пока просто картинки, а к формулам попозже вернусь.
…
(финал)
О, код на питоне! Написан вроде неплохо. Значит, теперь в формулы не обязательно вникать, можно же и отладчиком прогнать, если вдруг что. Облегчение.
(вместо заключения)
Многие авторы математических алгоритмов не прикладывают к своим статьям никакого кода, считая, что для программиста его математические выкладки должны быть очевидными. К сожалению, это редко когда так. Для программиста конструкция for(){for(){...}} более понятна, чем ? ?, даже если знать, что ? — это сумма, и несмотря на то, что она длиннее. Поэтому автору отдельный респект за код.
pchelintsev_an
09.12.2016 16:14Вот это хороший материал! Идеально подходит как для исследований, так и для лекций. А есть Интернет-ссылки на какие-нибудь работы автора в научных журналах, чтобы потом сослаться?
The_Freeman
09.12.2016 17:14Как я писал выше, Вы можете сослаться на книгу «Curve and Surface fitting with splines» автора Paul Dierckx. Я лишь модифицировал взятый оттуда алгоритм.
А лично у меня лишь одна статья в англоязычном научном журнале, и ссылаться на неё Вам вряд ли придётся =)
gudvinr
У вас есть возможность собрать пакет для PyPI?
hombit
Если я правильно понимаю, все уже сделано в
scipy
: https://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.UnivariateSpline.html