Интегрирование кусочно-линейной функции на неравномерной сет

[zcorvid]

Есть массив отсчётов, последовательность вещественных чисел $x_0, x_1, \dots , x_n$, ему соответствует массив времён $t_0, t_1, \dots , t_n$, нужно эти отсчёты как-то сгладить. Один из способов состоит в следующем: нарисуем пары точек $(t_i, x_i)$ на координатной плоскости $(t, x)$, рассмотрим их как как точки на графике некоторой (неизвестной нам) функции $x = x(t)$, соединим эти точки последовательно ломаной (то есть просто берётся кусочно линейная интерполяция). за пределы отрезка $[t_0, t_n]$ продолжим функцию константами: влево константой $x_0$, вправо $x_n$. Получившуюся кусочно-линейную функцию обозначим через $f(t)$. Получившийся график не очень красив, хочется его сгладить, для этого зададимся некоторым числом $\delta$, и определим функцию:

$$ F(t) = \int\limits_{t - \delta}^{t + \delta} f(t) dt $$

Задача состоит в том, чтобы (эффективно) посчитать массив значений $F(t_0), F(t_1), \dots , F(t_n)$.

Наивный алгоритм имеет асимптотику $O(n^2)$, если немного улучшить его бинарным поиском, можно добиться асимптотики $O(nlogn)$.

Набросок решения за $O(n)$.

Положим $\Delta t_i = t_i - t_{i - 1}$.

Рассмотрим интеграл кусочно-линейной функции по некоторому отрезку $[a, b]$, интеграл это площадь под графиком, значит это сумма площадей трапеций. Он распадается на две части, первая - между двумя точками из сетки $\{t_n\}$, и два маленьких кусочка по краям (первая часть это часть между крайними точками сетки, которые попадают в отрезок).

Интеграл по двум маленьким кусочкам считается за $O(1)$, а интеграл между узлами сетки можно подсчитать за $O(n)$, если предподсчитать кумулятивные суммы $S_k = \sum\limits_{i = 1}^{i < k} x_i (t_{i + 1} - t_{i - 1})$, тогда интеграл между узлами сетки $t_l, t_r$ выразится как $S_{l + 1} - S_{r} + x_l \Delta t_{l + 1} + x_r \Delta t_r$ (могу написать выкладку, как это получилось, подробно, если кратко - это получится, если в выражении для интеграла привести подобные члены по $x_i$).

Массив кумулятивных сумм считается за один проход, то есть $O(n)$, потом за второй проход считаем интеграл.

Пока не ясен один момент - как определять $l$ и $r$ для данных $a, b$, если делать это бинарным поиском, то итоговая асимптотика получится $O(nlogn)$. Скорее всего их можно определять адаптивно, используя информацию на предыдущих шагах, но я пока не придумал, как именно.

Вопрос - можно ли решить эту задачу за $O(n)$? И как её правильно решать? Если пытаться закодить вышеописанный алгоритм, получается какой-то монстр (походу я слишком криво его реализую). Возможно я изобретаю велосипед, и эта задача имеет решение каким-нибудь известным алгоритмом.

[Pphantom]

zcorvid в сообщении #1276036 писал(а):

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

Похоже. Посмотрите на аппроксимационные сплайны (например, кубические), при решении системы методом прогонки сложность получится $O(n)$.

[_Ivana]

Действительно, вы изобретаете целое семейство велосипедов. Уважаемый супермодератор предлагает глобальный сплайн, через решение линейной системы уравнений порядка , равного числу отсчетов, но через метод прогонки, который существенно снижает ужас вычислений. Я же со своей стороны хочу предложить любой из 100500 локальных сплайнов (любой степени гладкости, на выбор), которые для своего расчета на каждом интервале потребуют всего несколько сложений/умножений (при оптимизированном расчете). Начать можете хоть с самого простейшего - Катмулла-Рома (разновидности Эрмита) - расчет тривиален, а результат вас скорее всего устроит. Если нет - выбор локальных сплайнов весьма велик.

[Евгений Машеров]

Область изменения t распадается на два несвязных подмножества. Одно составляют значения t, удалённые от узлов на величину $2\delta$, второе $2\delta$-окрестности узлов. В первом подмножестве результат сглаживания кусочно-линейной функции будет совпадать с ней самой (предположу, что множитель $\frac 1 {2\delta}$ перед интегралом в выражении для F(t) случайно опущен, без него сглаживание константы её изменит). Во втором будет работать квадратичная интерполяция. Провести параболу можно по трём точкам. Например, по значению функции в точке $t_i-2\delta$, в точке $t_i+2\delta$ (обе и по значению требуемого интеграла в точке $t_i$ (который считается без интегрирования, просто зная наклоны линейных аппроксимаций слева и справа).
Отмечу, что предполагается, что $\delta<(t_i-t_{i+1})$ для всех i, в противном случае надо использовать обычные сплайны.

[zcorvid]

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

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

Множитель $\frac{1}{2\delta}$ действительно пропущен. А насчёт $\delta$ как раз наоборот, предполагается, что оно захватит некоторое количество соседних узлов. Идея в том, что можно было бы взять просто среднее 2 (4, 6, 8, ...) соседних отсчётов, и принять его за значение в точке, но при этом хочется, чтобы более удалённые по времени точки слабее влияли на значения (к примеру, если разница между точками более 10 секунд, то колебание между ними скорее всего обусловлено объективными причинами, а не ошибками измерения, и сглаживать это возмущение не следует, а если точки идут очень плотно с разницей в $0.1$ секунду, и быстро колеблются, то это уже скорее всего погрешность, так как реальный объект дёргаться не может, и в этом случае для усреднения можно взять побольше точек).

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

[_Ivana]

У меня ощущение, что вы описываете типичный цифровой ФНЧ. Разве что неравномерная сетка чуть разнообразит задачу.

[arseniiv]

Если дискретный анализ обобщён на функции на неравномерной сетке (а складывается впечатление, что давно уже написана по нему какая-нибудь монография) и там можно говорить о свёртке, можно было бы взять свёртку $f$ с чем-нибудь. (Упоминаемый _Ivana фильтр получится, если ядро свёртки зануляется слева от некоторой точки.)

[Евгений Машеров]

Да. Передискретизация+ФНЧ.

[Pphantom]

zcorvid в сообщении #1276386 писал(а):

Мне казалось, сплайны в точках сетки будут в точности равны заданным значениям $x_i$,

Это если они интерполяционные. А предлагаются аппроксимационные.

Ну или действительно ФНЧ.

[slavav]

Интеграл от кусочно-линейной функции есть кусочно-квадратичная функция. Строится за линейное время (заметанием, если использовать высокий слог). Значение этой функции в точке вычисляется за логарифм. Интеграл по отрезку есть разница двух значений интеграла на концах отрезка.
Сложность: $O(N)$ - предобработка, $O(\log N)$ - вычисление в точке. Можно показать что обе оценки оптимальные.

Полученную функцию вы табулируете. Наивная табуляция займёт время $O(K\log N)$, где $K$ - число точек табуляции. Так как последовательность точек табуляции возрастающая, то можно сэкономить на двоичном поиске. Точки отыскиваются последовательно линейным поиском (это тоже вариант заметания). Тогда сложность табуляции будет $O(\log N + K + M)$, где $M$ - число узлов изначальной функции на интервале табуляции, $M \leqslant N$. Если вы табулируете всю функцию, то $M = N$. В результате сложность будет $O(N + K)$.

-- 19.12.2017, 16:09 --

В вашем варианте ещё и $K = N$. Сложность будет $O(N)$.

[B@R5uk]

zcorvid в сообщении #1276386 писал(а):

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

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

[Евгений Машеров]

Я буду первый, кто произнесёт "Фильтр Калмана"?