Борьба с мельницами — 1: интерполяционные сплайны

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

— Где вы видите великанов? — спросил Санчо Панса.
— Да вон они, с громадными руками, — отвечал его господин. — У некоторых из них длина рук достигает почти двух миль.
— Помилуйте, сеньор, — возразил Санчо, — то, что там виднеется, вовсе не великаны, а ветряные мельницы; то же, что вы принимаете за их руки, — это крылья: они кружатся от ветра и приводят в движение мельничные жернова.
— Сейчас видно неопытного искателя приключений, — заметил Дон Кихот, — это великаны. И если ты боишься, то отъезжай в сторону и помолись, а я тем временем вступлю с ними в жестокий и неравный бой…

Однажды один мой друг попросил меня промоделировать ему простой и малозатратный вариант ресемплинга (изменения частоты дискретизации) цифрового сигнала. Особых требований к линейности АЧХ и ФЧХ не было, а желание экономить время и память было, поэтому часть методов типа полифазных интерполирующих банков FIR-фильтров сразу была исключена из рассмотрения, а основное внимание уделено методам локальной полиномиальной интерполяции. Полиномы быстро и легко считать, для хранения их коэффициентов требуется мало памяти, и с помощью одного набора коэффициентов можно посчитать любые значения внутри исходного интервала дискретизации.

Знакомимся с противником — реализация Фарроу

Известно, что через любые N точек (отсюда и в дальнейшем для определенности мы будем рассматривать равномерную сетку — с одинаковым шагом, хотя для многих описываемых техник и методов это не является обязательным условием), так вот, через любые N точек с разными координатами абсцисс можно провести единственный полином N-1 порядка, однозначно определяемый набором N коэффициентов (почитал, что LaTex на хабре не напрямую, а через картинки, которые тяжелы и недолговечны, поэтому буду обходиться только верхними-нижними индексами в html-тегах):

f(t) = a_n tⁿ + a_n-1 t^n-1 +…..+ a₁ t + a₀

Для удобства расчетов нормируем нашу равномерную сетку аргумента к целочисленным значениям, на примере четырех точек выберем удобный интервал аргумента (Фарроу выбрал от [-2, 1] с центральным интервалом от [-1, 0]) и запишем систему линейных уравнений из условия прохождения полинома через наши 4 точки:

f(i) = y_i
где i = {-2, -1, 0, 1}

Решать можно как угодно, хоть Гауссом матрицу обращать, но Фарроу оптимизировал количество операций для расчета коэффициентов этого полинома:

a₃ = (y₁ — y_-2)/6 + (y_-1 — y₀)/2
a₁ = (y₁ — y_-1)/2 — a₃
a₂ = y₁ — y₀ — a₁ — a₃
a₀ = y₀

Итого операций:

• 1 умножение на константу (1/6)
• 2 деления на 2 (считается быстро и в целочисленной арифметике и в формате с плавающей точкой)
• 8 сложений/вычитаний

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

Первый претендент — сплайн Катмулла-Рома

Можно заметить, что как сам искомый полином, так и все его производные являются линейными комбинациями неизвестных коэффициентов полинома. Поэтому мы можем задавать условия равенства производных любого порядка в узлах сетки заданным значениям, и получать линейную систему уравнений на коэффициенты полинома. Локальный кубический сплайн Эрмита строится по четырем условиям: значениям полинома и его первой производной в краях интервала. А получить приближенные значения этих производных в нужных точках мы можем, продифференцировав интерполирующий полином Лагранжа, проходящий через них. Например, нормируем аргументы четырех точек к интервалу [-1, 2], возьмем три точки с аргументами {-1, 0, 1} проведем через них параболу, и рассчитаем ее коэффициенты. А поскольку нам нужно только значение первой производной этой параболы в центральной точке 0, то достаточно рассчитать ее коэффициент при первой степени аргумента, он и будет использоваться в качестве производной для построения сплайна. Производная в правой границе интервала рассчитывается аналогично. В результате нехитрого расчета получаем следующую систему уравнений для коэффициентов сплайна:

f(i) = y_i
f ‘(i) = (y_i+1 — y_i-1)/2
где i = {0, 1}

решение которой можно получить в следующем виде:

y ‘₀ = (y₁ — y_-1)/2
y ‘₁ = (y₂ — y₀)/2

a₀ = y₀
a₁ = y ‘₁
a₃ = y ‘₀ + y ‘₁ + 2(y₀ — y₁)
a₂ = y₁ — a₃ — a₁ — a₀

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

• 2 умножения/деления на 2
• 7 сложений/вычитаний

Этот частный случай сплайна Эрмита и носит название Катмулла-Рома.

Второй претендент — сплайн Эрмита по 6 точкам

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

y ‘₀ = ( (y₁ — y_-1) — (y₂ — y_-2) /8 ) 2/3
y ‘₁ = ( (y₂ — y₀ ) — (y₃ — y_-1) /8 ) 2/3

a₀ = y₀
a₁ = y ‘₁
a₃ = y ‘₀ + y ‘₁ + 2(y₀ — y₁)
a₂ = y₁ — a₃ — a₁ — a₀

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

• 1 умножение на константу (2/3)
• 2 умножения/деления на степень 2
• 9 сложений/вычитаний

Сравнение вариантов

По количеству операций сплайн Катмулла-Рома существенно выигрывает у Фарроу, Эрмит по 6 точкам чуть проигрывает (на одно сложение). Также стоит отметить, что оба сплайна имеют первый порядок гладкости (непрерывные нулевую и первую производные) вследствие своего построения, в отличие от Фарроу, который не обеспечивает непрерывность производной. Но остается вопрос точности аппроксимации функций этими сплайнами. В качестве тестовой функции был выбран синус, произведена его аппроксимация вышеприведенными вариантами сплайнов для различного шага сетки, который нормируется в относительную величину количества точек на период. Измерялась величина максимального отклонения от аппроксимируемой функции. Результаты представлены на графике, при выборе логарифмического масштаба по обеим осям зависимости получаются почти линейные, сходящиеся в одну точку при уменьшении частоты дискретизации и приближению ее к частоте Найквиста (2 точки на период аппроксимируемого тестового синуса).

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

Заклинание

Мы можем развить общую идею сплайна Эрмита на производные любых порядков. Например, построим сплайн по следующим условиям: его значения в краях интервала совпадают со значениями отсчетов в этих точках, и вторые производные в этих точках совпадают с их оценками, произведенными по параболам через тройки точек — то есть, все как у Катмулла-Рома, только вместо первых будем брать вторые производные — парабола, как полином второго порядка, вполне позволит нам их рассчитать. И тут мы с удивлением обнаружим, что построенный таким образом сплайн полностью совпадает с нашим уважаемым Фарроу, который на самом деле полином Лагранжа. Желающие могут проверить это, я даже напишу формулу оценки второй производной параболой по 3 точкам (хотя и это тривиальное упражнение):

y » _i = y_i+1 + y_i-1 — 2 y_i

Удачное совпадение? Случайно, потому что мало точек и мал порядок полинома? Численно проверив этот факт на разном количестве точек и порядке полиномов, я убедился в том, что заклинание работает:

Пусть мы имеем интерполяционный полином, построенный по некоторым значениям на четном количестве точек равномерной сетки

x_-k, x_-(k-1),…x₀, x₁,…x_k-1, x_k, x_k+1

Тогда все четные производные этого полинома в точке x₀ совпадают с производными тех же порядков интерполяционного полинома, построенного по точкам

x_-k, x_-(k-1),…x₀, x₁,…x_k-1, x_k

а в точке x₁ — полинома, построенного по точкам

x_-(k-1),…x₀, x₁,…x_k-1, x_k, x_k+1

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

x_-k, x_-(k-1),…0, x₁,…x_k-1, x_k

Существует ее единственный интерполяционный полином P₀(x). Добавим к этому набору точек еще одну x_k+1 (для доказательства не важно, добавляем мы ее справа, слева или вообще в середину и не в узлы исходной сетки). Интерполяционный полином на расширенном наборе точек при условии равномерной сетки можно представить в форме Ньютона

P₁(x) = P₀(x) + A(x-x_-k)(x-x_-(k-1))….(x-0)….(x-x_k-1)(x-x_k)

где A — некая константа. Действительно, от этой аддитивной добавки требуется нулевые значения в точках исходного полинома, а константа A определяется из условия равенства P₁(x_k+1) = y_k+1 — требуемое значение полинома в добавленной точке. В силу равномерности сетки -x_-i = x_i, значит

P₁(x) = P₀(x) + A(x²-x_k²)(x²-x_k-1²)….(x²-x₁²)(x-0)

Очевидно, что в правой части стоит сумма P₀(x) и полинома с ненулевыми коэффициентами только при нечетных степенях аргумента. Значит, коэффициенты при четных степенях у полиномов P₁(x) и P₀(x) совпадают, а все четные производные любого полинома в нуле определяются только соответствующим его коэффициентом при четной степени аргумента. Таким образом, все четные производные P₁(x) и P₀(x) совпадают.
Вышеприведенные рассуждения справедливы для любого расположения добавляемой точки x_k+1, в частности при добавлении ее справа или слева по исходной сетке. Таким образом, производные четных порядков в нуле у исходного полинома и полинома по расширенному набору точек совпадают. Ч.т.д.

Последний раунд — Лагранж по Фарроу с Лагранжем через производные

Поскольку интерполяционный полином минимального порядка единственен, то мы можем с учетом только что доказанного заклинания строить его через условия на значения четных производных в краях центрального интервала — например, полином 5-го порядка через 6 точек по 6 условиям: значения отсчетов и вторых и четвертых производных, рассчитанных по каждым пяти точкам. При этом пользоваться той же вышеописанной оптимизацией — расчет производных только в правой границе каждого интервала. В первоначальной постановке прохождения полинома через все точки эти инварианты были не очевидны. Применим это к нашему исходному Лагранжу по 4 точкам и оптимизируем количество операций. Приведу сразу код Matlab своего варианта, и варианта участника с ником _Anatoliy с форума электроникса, который независимо и самостоятельно побеждал Фарроу еще до моих экспериментов:

clf reset N = 10; x = 1:N; y = rand(1, N); h = 0.01; g = y(3) - y(2); d = g - y(2) + y(1); plot(x, y, 'or', 'LineWidth', 2); hold on; grid on; axis on; title(['Локальная интерполяция случайного набора точек', ... ' полиномом Лагранжа 3 степени, 2 алгоритма расчета.']); plot(x(1), y(1), 'b-', x(1), y(1), 'g:'); legend('точки', 'мой алгоритм', 'алгоритм Anatoliy'); for k = 2:(N-2)     % мой алгоритм: 1 умножение, 1 сдвиг, 5 сложений     c = d; e = g; g = y(k+2) - y(k+1); d = g - e;     b2 = c/2; b3 = (d - c)/6; b1 = e - b2 - b3;     t = 0:h:1; f = b3.*t.^3 + b2.*t.^2 + b1.*t + y(k);     plot(t+x(k), f, 'b-')          % алгоритм _Anatoliy: 1 умножение, 2 сдвига, 6 сложений            p = y(k) - y(k+1); q = (y(k+2) - y(k))/2;     a3 = ( (y(k+2) - y(k-1))/3 + p )/2;     a2 = p + q; a1 = q - a3;          t = -1:h:0; f = a3.*t.^3 + a2.*t.^2 + a1.*t + y(k+1);     plot(t+1+x(k), f, 'g:') end 

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

ссылка на оригинал статьи https://habrahabr.ru/post/282441/