В данной статье лирический герой бросает вызов оптимальной реализации классического полиномиального интерполятора Лагранжа (Фарроу), в процессе битвы случайно открывает и доказывает тривиальное никому не нужное математгическое заклинание, с помощью которого пытается потеснить противника, но по результатам всех раундов боя решением судей фиксируется ничья.
Однажды один мой друг попросил меня промоделировать ему простой и малозатратный вариант ресемплинга (изменения частоты дискретизации) цифрового сигнала. Особых требований к линейности АЧХ и ФЧХ не было, а желание экономить время и память было, поэтому часть методов типа полифазных интерполирующих банков FIR-фильтров сразу была исключена из рассмотрения, а основное внимание уделено методам локальной полиномиальной интерполяции. Полиномы быстро и легко считать, для хранения их коэффициентов требуется мало памяти, и с помощью одного набора коэффициентов можно посчитать любые значения внутри исходного интервала дискретизации.
Известно, что через любые N точек (отсюда и в дальнейшем для определенности мы будем рассматривать равномерную сетку — с одинаковым шагом, хотя для многих описываемых техник и методов это не является обязательным условием), так вот, через любые N точек с разными координатами абсцисс можно провести единственный полином N-1 порядка, однозначно определяемый набором N коэффициентов (почитал, что LaTex на хабре не напрямую, а через картинки, которые тяжелы и недолговечны, поэтому буду обходиться только верхними-нижними индексами в html-тегах):
f(t) = an tn + an-1 tn-1 +.....+ a1 t + a0
Для удобства расчетов нормируем нашу равномерную сетку аргумента к целочисленным значениям, на примере четырех точек выберем удобный интервал аргумента (Фарроу выбрал от [-2, 1] с центральным интервалом от [-1, 0]) и запишем систему линейных уравнений из условия прохождения полинома через наши 4 точки:
f(i) = yi
где i = {-2, -1, 0, 1}
Решать можно как угодно, хоть Гауссом матрицу обращать, но Фарроу оптимизировал количество операций для расчета коэффициентов этого полинома:
a3 = (y1 — y-2)/6 + (y-1 — y0)/2
a1 = (y1 — y-1)/2 — a3
a2 = y1 — y0 — a1 — a3
a0 = y0
Итого операций:
Данный полином используется для расчета значений, попадающих в центральный интервал исходных четырех точек: [-1, 0] (с соответствующим нормированием аргумента), для других интервалов аналогичным образом рассчитываются полиномы по окружающим их четырем точкам. Очень популярный интерполятор в узких кругах, заслуженный и почитаемый. И мне почему-то сразу захотелось придумать что-то, что может превзойти эту реализацию в каком-то смысле, с учетом количества операций для расчета.
Можно заметить, что как сам искомый полином, так и все его производные являются линейными комбинациями неизвестных коэффициентов полинома. Поэтому мы можем задавать условия равенства производных любого порядка в узлах сетки заданным значениям, и получать линейную систему уравнений на коэффициенты полинома. Локальный кубический сплайн Эрмита строится по четырем условиям: значениям полинома и его первой производной в краях интервала. А получить приближенные значения этих производных в нужных точках мы можем, продифференцировав интерполирующий полином Лагранжа, проходящий через них. Например, нормируем аргументы четырех точек к интервалу [-1, 2], возьмем три точки с аргументами {-1, 0, 1} проведем через них параболу, и рассчитаем ее коэффициенты. А поскольку нам нужно только значение первой производной этой параболы в центральной точке 0, то достаточно рассчитать ее коэффициент при первой степени аргумента, он и будет использоваться в качестве производной для построения сплайна. Производная в правой границе интервала рассчитывается аналогично. В результате нехитрого расчета получаем следующую систему уравнений для коэффициентов сплайна:
f(i) = yi
f '(i) = (yi+1 — yi-1)/2
где i = {0, 1}
решение которой можно получить в следующем виде:
y '0 = (y1 — y-1)/2
y '1 = (y2 — y0)/2
a0 = y0
a1 = y '1
a3 = y '0 + y '1 + 2(y0 — y1)
a2 = y1 — a3 — a1 — a0
Можем заметить, что при последовательной (потоковой) обработке входных данных производная левой границы интервала та же самая, что и производная правой границы предыдущего интервала, и значит на каждом интервале нам необходимо рассчитывать только одну производную — в правой границе интервала, а значение в левой брать из запомненной при предыдущем расчете в правой. С помощью такой экономии на спичках получаем следующее количество операций для расчета коэффициентов полинома:
Этот частный случай сплайна Эрмита и носит название Катмулла-Рома.
Если мы будем оценивать производные в краях центрального интервала не по трем, а по пяти точкам, то их расчет немного усложнится, но остальные формулы не поменяются:
y '0 = ( (y1 — y-1) — (y2 — y-2) /8 ) 2/3
y '1 = ( (y2 — y0 ) — (y3 — y-1) /8 ) 2/3
a0 = y0
a1 = y '1
a3 = y '0 + y '1 + 2(y0 — y1)
a2 = y1 — a3 — a1 — a0
С учетом предыдущего замечания про необходимость расчета только одной производной на каждом интервале, количество операций:
По количеству операций сплайн Катмулла-Рома существенно выигрывает у Фарроу, Эрмит по 6 точкам чуть проигрывает (на одно сложение). Также стоит отметить, что оба сплайна имеют первый порядок гладкости (непрерывные нулевую и первую производные) вследствие своего построения, в отличие от Фарроу, который не обеспечивает непрерывность производной. Но остается вопрос точности аппроксимации функций этими сплайнами. В качестве тестовой функции был выбран синус, произведена его аппроксимация вышеприведенными вариантами сплайнов для различного шага сетки, который нормируется в относительную величину количества точек на период. Измерялась величина максимального отклонения от аппроксимируемой функции. Результаты представлены на графике, при выборе логарифмического масштаба по обеим осям зависимости получаются почти линейные, сходящиеся в одну точку при уменьшении частоты дискретизации и приближению ее к частоте Найквиста (2 точки на период аппроксимируемого тестового синуса).
На графике видно, что сплайн Катмулла-Рома уступает Фарроу по точности, а последний уступает Эрмиту по 6 точкам, при практически одинаковом количестве требуемых для расчета операций. Одержать чистую победу не удалось, но в принципе результат имхо неплохой.
Мы можем развить общую идею сплайна Эрмита на производные любых порядков. Например, построим сплайн по следующим условиям: его значения в краях интервала совпадают со значениями отсчетов в этих точках, и вторые производные в этих точках совпадают с их оценками, произведенными по параболам через тройки точек — то есть, все как у Катмулла-Рома, только вместо первых будем брать вторые производные — парабола, как полином второго порядка, вполне позволит нам их рассчитать. И тут мы с удивлением обнаружим, что построенный таким образом сплайн полностью совпадает с нашим уважаемым Фарроу, который на самом деле полином Лагранжа. Желающие могут проверить это, я даже напишу формулу оценки второй производной параболой по 3 точкам (хотя и это тривиальное упражнение):
y '' i = yi+1 + yi-1 — 2 yi
Удачное совпадение? Случайно, потому что мало точек и мал порядок полинома? Численно проверив этот факт на разном количестве точек и порядке полиномов, я убедился в том, что заклинание работает:
Пусть мы имеем интерполяционный полином, построенный по некоторым значениям на четном количестве точек равномерной сетки
x-k, x-(k-1),...x0, x1,...xk-1, xk, xk+1
Тогда все четные производные этого полинома в точке x0 совпадают с производными тех же порядков интерполяционного полинома, построенного по точкам
x-k, x-(k-1),...x0, x1,...xk-1, xk
а в точке x1 — полинома, построенного по точкам
x-(k-1),...x0, x1,...xk-1, xk, xk+1
Доказательство. Рассмотрим функцию на равномерной сетке из нечетного количества точек — сделаем смещение аргумента, приводящее аргумент центрального узла в 0
x-k, x-(k-1),...0, x1,...xk-1, xk
Существует ее единственный интерполяционный полином P0(x). Добавим к этому набору точек еще одну xk+1 (для доказательства не важно, добавляем мы ее справа, слева или вообще в середину и не в узлы исходной сетки). Интерполяционный полином на расширенном наборе точек при условии равномерной сетки можно представить в форме Ньютона
P1(x) = P0(x) + A(x-x-k)(x-x-(k-1))....(x-0)....(x-xk-1)(x-xk)
где A — некая константа. Действительно, от этой аддитивной добавки требуется нулевые значения в точках исходного полинома, а константа A определяется из условия равенства P1(xk+1) = yk+1 — требуемое значение полинома в добавленной точке. В силу равномерности сетки -x-i = xi, значит
P1(x) = P0(x) + A(x2-xk2)(x2-xk-12)....(x2-x12)(x-0)
Очевидно, что в правой части стоит сумма P0(x) и полинома с ненулевыми коэффициентами только при нечетных степенях аргумента. Значит, коэффициенты при четных степенях у полиномов P1(x) и P0(x) совпадают, а все четные производные любого полинома в нуле определяются только соответствующим его коэффициентом при четной степени аргумента. Таким образом, все четные производные P1(x) и P0(x) совпадают.
Вышеприведенные рассуждения справедливы для любого расположения добавляемой точки xk+1, в частности при добавлении ее справа или слева по исходной сетке. Таким образом, производные четных порядков в нуле у исходного полинома и полинома по расширенному набору точек совпадают. Ч.т.д.
Поскольку интерполяционный полином минимального порядка единственен, то мы можем с учетом только что доказанного заклинания строить его через условия на значения четных производных в краях центрального интервала — например, полином 5-го порядка через 6 точек по 6 условиям: значения отсчетов и вторых и четвертых производных, рассчитанных по каждым пяти точкам. При этом пользоваться той же вышеописанной оптимизацией — расчет производных только в правой границе каждого интервала. В первоначальной постановке прохождения полинома через все точки эти инварианты были не очевидны. Применим это к нашему исходному Лагранжу по 4 точкам и оптимизируем количество операций. Приведу сразу код Matlab своего варианта, и варианта участника с ником _Anatoliy с форума электроникса, который независимо и самостоятельно побеждал Фарроу еще до моих экспериментов:
Как видно, оба кода содержат меньшее количество операций для расчета коэффициентов полинома, так что в плане спорта я считаю, что мы оба одержали победы. И хотя высказываются мнения, что это экономия на спичках, что в наше время, когдакосмические корабли бороздят везде есть математические сопроцессоры с интегрированными операциями с числами с плавающей точкой, скорость операций сравнима со скоростью перемещения содержимого регистров и тем более чтения/записи памяти для сохранения/восстановления рассчитанных значений — подобные оптимизации ничего не дают, у меня все равно есть некое чувство удовлетворенности результатами, несмотря на их слабую практическую значимость.
— Где вы видите великанов? — спросил Санчо Панса.
— Да вон они, с громадными руками, — отвечал его господин. — У некоторых из них длина рук достигает почти двух миль.
— Помилуйте, сеньор, — возразил Санчо, — то, что там виднеется, вовсе не великаны, а ветряные мельницы; то же, что вы принимаете за их руки, — это крылья: они кружатся от ветра и приводят в движение мельничные жернова.
— Сейчас видно неопытного искателя приключений, — заметил Дон Кихот, — это великаны. И если ты боишься, то отъезжай в сторону и помолись, а я тем временем вступлю с ними в жестокий и неравный бой…
Однажды один мой друг попросил меня промоделировать ему простой и малозатратный вариант ресемплинга (изменения частоты дискретизации) цифрового сигнала. Особых требований к линейности АЧХ и ФЧХ не было, а желание экономить время и память было, поэтому часть методов типа полифазных интерполирующих банков FIR-фильтров сразу была исключена из рассмотрения, а основное внимание уделено методам локальной полиномиальной интерполяции. Полиномы быстро и легко считать, для хранения их коэффициентов требуется мало памяти, и с помощью одного набора коэффициентов можно посчитать любые значения внутри исходного интервала дискретизации.
Знакомимся с противником — реализация Фарроу
Известно, что через любые N точек (отсюда и в дальнейшем для определенности мы будем рассматривать равномерную сетку — с одинаковым шагом, хотя для многих описываемых техник и методов это не является обязательным условием), так вот, через любые N точек с разными координатами абсцисс можно провести единственный полином N-1 порядка, однозначно определяемый набором N коэффициентов (почитал, что LaTex на хабре не напрямую, а через картинки, которые тяжелы и недолговечны, поэтому буду обходиться только верхними-нижними индексами в html-тегах):
f(t) = an tn + an-1 tn-1 +.....+ a1 t + a0
Для удобства расчетов нормируем нашу равномерную сетку аргумента к целочисленным значениям, на примере четырех точек выберем удобный интервал аргумента (Фарроу выбрал от [-2, 1] с центральным интервалом от [-1, 0]) и запишем систему линейных уравнений из условия прохождения полинома через наши 4 точки:
f(i) = yi
где i = {-2, -1, 0, 1}
Решать можно как угодно, хоть Гауссом матрицу обращать, но Фарроу оптимизировал количество операций для расчета коэффициентов этого полинома:
a3 = (y1 — y-2)/6 + (y-1 — y0)/2
a1 = (y1 — y-1)/2 — a3
a2 = y1 — y0 — a1 — a3
a0 = y0
Итого операций:
- 1 умножение на константу (1/6)
- 2 деления на 2 (считается быстро и в целочисленной арифметике и в формате с плавающей точкой)
- 8 сложений/вычитаний
Данный полином используется для расчета значений, попадающих в центральный интервал исходных четырех точек: [-1, 0] (с соответствующим нормированием аргумента), для других интервалов аналогичным образом рассчитываются полиномы по окружающим их четырем точкам. Очень популярный интерполятор в узких кругах, заслуженный и почитаемый. И мне почему-то сразу захотелось придумать что-то, что может превзойти эту реализацию в каком-то смысле, с учетом количества операций для расчета.
Первый претендент — сплайн Катмулла-Рома
Можно заметить, что как сам искомый полином, так и все его производные являются линейными комбинациями неизвестных коэффициентов полинома. Поэтому мы можем задавать условия равенства производных любого порядка в узлах сетки заданным значениям, и получать линейную систему уравнений на коэффициенты полинома. Локальный кубический сплайн Эрмита строится по четырем условиям: значениям полинома и его первой производной в краях интервала. А получить приближенные значения этих производных в нужных точках мы можем, продифференцировав интерполирующий полином Лагранжа, проходящий через них. Например, нормируем аргументы четырех точек к интервалу [-1, 2], возьмем три точки с аргументами {-1, 0, 1} проведем через них параболу, и рассчитаем ее коэффициенты. А поскольку нам нужно только значение первой производной этой параболы в центральной точке 0, то достаточно рассчитать ее коэффициент при первой степени аргумента, он и будет использоваться в качестве производной для построения сплайна. Производная в правой границе интервала рассчитывается аналогично. В результате нехитрого расчета получаем следующую систему уравнений для коэффициентов сплайна:
f(i) = yi
f '(i) = (yi+1 — yi-1)/2
где i = {0, 1}
решение которой можно получить в следующем виде:
y '0 = (y1 — y-1)/2
y '1 = (y2 — y0)/2
a0 = y0
a1 = y '1
a3 = y '0 + y '1 + 2(y0 — y1)
a2 = y1 — a3 — a1 — a0
Можем заметить, что при последовательной (потоковой) обработке входных данных производная левой границы интервала та же самая, что и производная правой границы предыдущего интервала, и значит на каждом интервале нам необходимо рассчитывать только одну производную — в правой границе интервала, а значение в левой брать из запомненной при предыдущем расчете в правой. С помощью такой экономии на спичках получаем следующее количество операций для расчета коэффициентов полинома:
- 2 умножения/деления на 2
- 7 сложений/вычитаний
Этот частный случай сплайна Эрмита и носит название Катмулла-Рома.
Второй претендент — сплайн Эрмита по 6 точкам
Если мы будем оценивать производные в краях центрального интервала не по трем, а по пяти точкам, то их расчет немного усложнится, но остальные формулы не поменяются:
y '0 = ( (y1 — y-1) — (y2 — y-2) /8 ) 2/3
y '1 = ( (y2 — y0 ) — (y3 — y-1) /8 ) 2/3
a0 = y0
a1 = y '1
a3 = y '0 + y '1 + 2(y0 — y1)
a2 = y1 — a3 — a1 — a0
С учетом предыдущего замечания про необходимость расчета только одной производной на каждом интервале, количество операций:
- 1 умножение на константу (2/3)
- 2 умножения/деления на степень 2
- 9 сложений/вычитаний
Сравнение вариантов
По количеству операций сплайн Катмулла-Рома существенно выигрывает у Фарроу, Эрмит по 6 точкам чуть проигрывает (на одно сложение). Также стоит отметить, что оба сплайна имеют первый порядок гладкости (непрерывные нулевую и первую производные) вследствие своего построения, в отличие от Фарроу, который не обеспечивает непрерывность производной. Но остается вопрос точности аппроксимации функций этими сплайнами. В качестве тестовой функции был выбран синус, произведена его аппроксимация вышеприведенными вариантами сплайнов для различного шага сетки, который нормируется в относительную величину количества точек на период. Измерялась величина максимального отклонения от аппроксимируемой функции. Результаты представлены на графике, при выборе логарифмического масштаба по обеим осям зависимости получаются почти линейные, сходящиеся в одну точку при уменьшении частоты дискретизации и приближению ее к частоте Найквиста (2 точки на период аппроксимируемого тестового синуса).
На графике видно, что сплайн Катмулла-Рома уступает Фарроу по точности, а последний уступает Эрмиту по 6 точкам, при практически одинаковом количестве требуемых для расчета операций. Одержать чистую победу не удалось, но в принципе результат имхо неплохой.
Заклинание
Мы можем развить общую идею сплайна Эрмита на производные любых порядков. Например, построим сплайн по следующим условиям: его значения в краях интервала совпадают со значениями отсчетов в этих точках, и вторые производные в этих точках совпадают с их оценками, произведенными по параболам через тройки точек — то есть, все как у Катмулла-Рома, только вместо первых будем брать вторые производные — парабола, как полином второго порядка, вполне позволит нам их рассчитать. И тут мы с удивлением обнаружим, что построенный таким образом сплайн полностью совпадает с нашим уважаемым Фарроу, который на самом деле полином Лагранжа. Желающие могут проверить это, я даже напишу формулу оценки второй производной параболой по 3 точкам (хотя и это тривиальное упражнение):
y '' i = yi+1 + yi-1 — 2 yi
Удачное совпадение? Случайно, потому что мало точек и мал порядок полинома? Численно проверив этот факт на разном количестве точек и порядке полиномов, я убедился в том, что заклинание работает:
Пусть мы имеем интерполяционный полином, построенный по некоторым значениям на четном количестве точек равномерной сетки
x-k, x-(k-1),...x0, x1,...xk-1, xk, xk+1
Тогда все четные производные этого полинома в точке x0 совпадают с производными тех же порядков интерполяционного полинома, построенного по точкам
x-k, x-(k-1),...x0, x1,...xk-1, xk
а в точке x1 — полинома, построенного по точкам
x-(k-1),...x0, x1,...xk-1, xk, xk+1
Доказательство. Рассмотрим функцию на равномерной сетке из нечетного количества точек — сделаем смещение аргумента, приводящее аргумент центрального узла в 0
x-k, x-(k-1),...0, x1,...xk-1, xk
Существует ее единственный интерполяционный полином P0(x). Добавим к этому набору точек еще одну xk+1 (для доказательства не важно, добавляем мы ее справа, слева или вообще в середину и не в узлы исходной сетки). Интерполяционный полином на расширенном наборе точек при условии равномерной сетки можно представить в форме Ньютона
P1(x) = P0(x) + A(x-x-k)(x-x-(k-1))....(x-0)....(x-xk-1)(x-xk)
где A — некая константа. Действительно, от этой аддитивной добавки требуется нулевые значения в точках исходного полинома, а константа A определяется из условия равенства P1(xk+1) = yk+1 — требуемое значение полинома в добавленной точке. В силу равномерности сетки -x-i = xi, значит
P1(x) = P0(x) + A(x2-xk2)(x2-xk-12)....(x2-x12)(x-0)
Очевидно, что в правой части стоит сумма P0(x) и полинома с ненулевыми коэффициентами только при нечетных степенях аргумента. Значит, коэффициенты при четных степенях у полиномов P1(x) и P0(x) совпадают, а все четные производные любого полинома в нуле определяются только соответствующим его коэффициентом при четной степени аргумента. Таким образом, все четные производные P1(x) и P0(x) совпадают.
Вышеприведенные рассуждения справедливы для любого расположения добавляемой точки xk+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
Как видно, оба кода содержат меньшее количество операций для расчета коэффициентов полинома, так что в плане спорта я считаю, что мы оба одержали победы. И хотя высказываются мнения, что это экономия на спичках, что в наше время, когда
saproman88
Сейчас нет возможности вечером найду рисунки.
Посмотрите в монографии: В.Ф. Снигирев. Применение сплайнов для задания обводов летательных аппаратов.
Там приводится пример интересной функции для тестирования сплайнов.
По моему она еще есть в книге: «Конструкция самолетов» авторы О.А. Гребеньков, В.П. Гоголин, А.И. Осокин, В.Ф. Снигирев, В.Г. Шатаев,
IIvana
Спасибо, попробую поискать ее в инете. Но может это спецлитература и не оцифрована. Может речь идет о функции Рунге, которая приводится во всех примерах, когда хотят показать осцилляции полиномов Лагранжа высоких порядков, и на основании этого делают вывод, что «полиномы высоких порядков непригодны для интерполяции», хотя это совсем не так.
saproman88
Завьялов, Квасов, Мирошниченко. Методы сплайн функций, стр. 193. Формула этой функции. а на 191 или 192 странице вид самой функции. Эта книга точно оцифрована и ее можно найти.
IIvana
А эта книга у меня есть, и даже неоднократно читана. Скорее всего о функции Рунге идет речь, вечером посмотрю. Но нет никаких проблем при ее интерполяции, если только не пытаться интерполировать ее одним единственным полиномом по многим точкам, а если использовать не глобальные, а локальные методы (все перечисленные в статье — локальны), тогда все хорошо.
IIvana
Посмотрел, там приводится пример функции подобной ступеньке Хэвисайда, с крутым фронтом спада, и интерполяция ее глобальными кубическим и рациональным сплайнами по 10 точкам, с демонстрацией преимущества последнего. Конечно, аргумент в знаменателе дает больше гибкости.