[Alexey Lukin 0]

Ну, следовало сказать, что интеграл несобственный.

[Xenia 42]

Лежандр рулит! :)

 

При длине вектора 15:

 

    1.0000   -0.0000    0.0714         0    0.0825         0
   -0.0000    0.3810         0    0.0777         0    0.0917
    0.0714         0    0.2628    0.0000    0.0870         0
         0    0.0777    0.0000    0.2187         0    0.0984
    0.0825         0    0.0870         0    0.2006    0.0000
         0    0.0917         0    0.0984    0.0000    0.1938

Шероховатости в виде внедиагональных элементов видны, но уже не такие большие, как у Чебышева, и могут быть списаны на грубость дискретизации (15 элементов это действительно мало).

Проверяем эту гипотезу, увеличивая число элементов в векторе до 101:

 

    1.0000   -0.0000    0.0100    0.0000    0.0102   -0.0000
   -0.0000    0.3400    0.0000    0.0101    0.0000    0.0104
    0.0100    0.0000    0.2081   -0.0000    0.0103   -0.0000
    0.0000    0.0101   -0.0000    0.1517   -0.0000    0.0106
    0.0102    0.0000    0.0103   -0.0000    0.1206    0.0000
   -0.0000    0.0104   -0.0000    0.0106    0.0000    0.1009

Увеличиваем число элементов в векторе до 1001:

 

    1.0000   -0.0000    0.0010   -0.0000    0.0010    0.0000
   -0.0000    0.3340   -0.0000    0.0010    0.0000    0.0010
    0.0010   -0.0000    0.2008    0.0000    0.0010   -0.0000
   -0.0000    0.0010    0.0000    0.1437   -0.0000    0.0010
    0.0010    0.0000    0.0010   -0.0000    0.1120    0.0000
    0.0000    0.0010   -0.0000    0.0010    0.0000    0.0918

Увеличиваем число элементов в векторе до 10001:

 

    1.0000   -0.0000    0.0001   -0.0000    0.0001   -0.0000
   -0.0000    0.3334   -0.0000    0.0001   -0.0000    0.0001
    0.0001   -0.0000    0.2001   -0.0000    0.0001    0.0000
   -0.0000    0.0001   -0.0000    0.1429    0.0000    0.0001
    0.0001   -0.0000    0.0001    0.0000    0.1112    0.0000
   -0.0000    0.0001    0.0000    0.0001    0.0000    0.0910

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

[_Ivana 0]

AndrewN а если при определенных парах полиномов с весом первообразная в +-1 будет в бесконечность улетать? Вы же написали интересные любопытные вещи, давайте не будем спорить об имхо несущественных деталях.

 

Xenia рад что у Вас все получилось да ещё сразу проверилось моделированием.

Изменено 24 июля, 2012 пользователем _Ivana

[andrewn 0]

Ну, следовало сказать, что интеграл несобственный

Первообразная-то все равно без особенностей по краям. Вычислили и вычли, а несобственность (пределы eps слева-справа) держим в уме.

 

 

 

Увеличиваем число элементов в векторе до ...

Такой метод можно назвать "переверну Дарбу в гробу"...

 

 

[_Ivana 0]

ТС нужна была система ортогональных полиномов в единичным весом - она их нашла. Для чего и как применять - это уже второй вопрос. Можно ли применить более оптимальные методы для решения этой неизвестной авторской задачи - уже третий. А "Увеличиваем число элементов в векторе до ..." в этом контексте - не более чем метод проверки базиса на ортогональность, который кстати работает.

[andrewn 0]

который кстати работает

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

 

 

[Xenia 42]

Такой метод можно назвать "переверну Дарбу в гробу"...

 

Отчего же? :) Интеграл изначально определялся как предел сумм при бесконечном увеличении разбиения отрезка на части. Это уж потом привыкли к разным аналитическим выражениям, черпаемых в справочниках :).

 

А я, увеличивая длину вектора, представляющего собой разбиение отрезка [-1,+1], делала буквально то самое, что определяется, как интеграл. Правда, до бесконечности я не дошла, а остановилась на 1/10000 от ширины интервала, но с моими 4-мя знаками после запятой этот предел вполне достаточен, поскольку большей точности при таком выводе результата я все равно не увижу.

 

На практике же мне нужны именно суммы, а не интегралы в их аналитическом выражении. Поскольку исходные данные, которые мне придется переводить в базис Лежандра, исключительно табличные (это данные от АЦП). И для этого перевода мне придется находить скалярные произведения векторов длиной gap. Если базисные вектора достаточно хорошо ортогональны, то можно было бы не раскручивать МНК с обращением матрицы, а просто поочередно перемножать вектор данных на вектора ортогонального базиса. Поэтому в данной задаче меня больше всего волнует именно та погрешность, которая возникнет, если стану работать с базисом, как ортогонльным, в условиях, когда дискретизация представления базиса в векторах реально дает ортогонализацию с ошибкой.

[andrewn 0]

А я, увеличивая длину вектора, представляющего собой разбиение отрезка [-1,+1], делала буквально то самое, что определяется, как интеграл. Правда, до бесконечности я не дошла, а остановилась на 1/10000 от ширины интервала, но с моими 4-мя знаками после запятой этот предел вполне достаточен, поскольку большей точности при таком выводе результата я все равно не увижу.

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

 

Ну в самом деле, на любом компьютере сумму гармонического ряда можно сосчитать как два пальца об асфальт, просто потому что а) суммирование элементарно приходится оборвать и б) начиная с некоторого n члены просто перестают давать вклад в бегущую сумму. "Но мы-то знаем", что ряд расходится!

 

 

 

исходные данные, [...] мне придется переводить в базис Лежандра

Арифметика мне понятна, вы считаете коэффициенты разложения по базису. Но что вас заставляет использовать такой неудобный базис? В этом какая-то физика замешана?

 

 

[Xenia 42]

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

Так я и не пыталась решать глобальные задачи типа сходимости неизвестных рядов. Суть была в том, что полиномы Чебышева не образуют ортогональную (в моем понимании!) систему, а в полиномы Эрмита ее образуют. За что еще раз спасибо _Ivana, что надоумил.

 

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

 

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

 

Арифметика мне понятна, вы считаете коэффициенты разложения по базису. Но что вас заставляет использовать такой неудобный базис? В этом какая-то физика замешана?

Почему вдруг этот базис неудобный? А какой тогда удобный? :)

 

Физически задача у меня довольно тривиальная - регрессия. Т.е. я хочу "впендюрить" отрезки экспериментальной кривой в базис невысокой размерности. Так очень часто поступают, аппроксимируя экспериментальные точки прямой, параболой, реже полиномами более высоких порядков. Я тоже так делала, но недовольна сильной корреляцией между собой степенных полиномов, из-за чего обращаемая матрица становится близкой к матрице Гилберта и обращается с ошибками.

 

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

 

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

 

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

[andrewn 0]

в базис невысокой размерности

Ага. А, кстати, можно ли посмотреть на пару-тройку файлов с измерениями - отправьте на е-мейл?

 

[Xenia 42]

А, кстати, можно ли посмотреть на пару-тройку файлов с измерениями - отправьте на е-мейл?

 

Посмотреть нельзя! :) Однако тут дело не в конкретных данных, а в самой "странно поставленной" задаче. Если хотите, то в общих чертах объясню ситуацию.

 

Ползет лента самописца... Только по нынешнему времени это не лента, а высокоразрядный АЦП данные снимает и на экране компьютера рисует. Ну и, как водится, сигнал шумит. Если бы это самописец был, то он бы толстую линию рисовал, т.к. перо у него ходуном бы ходило (из-за своей инерционности обычные бумажные самописцы на высокочастотный шум не реагируют, но АЦП его отлично видит).

 

Только не торопитесь ставить диагноз - "Так вам шум фильтровать надо? Мы это и без всяких полиномов смастырим!" :) Нет, не всё тут тривиально, а потому наберитесь терпения дослушать моё объяснение до конца. Так вот "задним числом" по тем сырым данным ползет окно и аппроксимирует точки/значения, в него попавшие, полиномом какой-то степени (математически это зовется регрессией). Ширина окна и степень полинома выбираются не с потолка, а исходя из определенных соображений, и в процессе продвижения окна могут автоматически изменяться в зависимости от того, насколько удачно проходит аппроксимация. Эти соображения я пока опущу. А дальше строго по статистическим критериям (там и t-критерий Стьдента и хи-квадрат пробабилити) оценивается, входит ли центральная точка в окне в доверительный интервал. Если входит, то подлежит замене на значение полинома в этой точке. Т.е. в этом случае имеет место типичное скользящее сглаживание по N-точкам полиномом n-ой степени. Но это лишь пока тишь да гладь, а шум укладывается в пределы статистической нормы. Но стоит какой-то точке или нескольким точкам не вписаться в "валютный коридор" :), то их аппроксимировать уже нельзя! И рисуются такие точки, как есть, без всякого сглаживания. Такие сигналы квалифицируются, как выбросы/феномены, и их изучать придут большие ученые :). И горе тому, кто сгладит выброс! :)

[andrewn 0]

Посмотреть нельзя

Ага. 3 сигмы или две?

А если вместо странных полиномов использовать косинусное преобразование для сигнала - модели то толковой для него всё рано нету...

И почему посмотреть-то нельзя - я по временному ряду кита от стаи селёдки все равно не отличу, а любопытство заело...

 

 

[Xenia 42]

Ага. 3 сигмы или две?

Две сигмы ("the Inverse Student-t Probability Distribution Function (two-tail) for 5% error").

 

А если вместо странных полиномов использовать косинусное преобразование для сигнала - модели то толковой для него всё рано нету...

Я уже на этот счет заикалась. Косинусное, как и обычное Фурье-преобразование плохо аппроксимирует импульные скачки. Например, если вдруг "дернулось" (короткий бросок напряжения), то нет такой частоты гармоники, которая бы удовлеворительно его аппроксимировала. А в результате произойдет следующее - получится целый спектр разночастотных гармоник, имеющих синхроннизацию "выпуклости" на этом месте. С тем, значит, расчетом, чтобы на месте импульса гармоники сложились в одной полярности, а по бокам друг дружку погасили. А если я ограничу размер базиса первыми n-гармониками, то после апроксимации получу по всему полю периодическую волну с часттой младшей гармоники базиса, которую некому будет погасить. Прошу прощения за слог, но лучше объяснить на словах это являение я не могу. А вот полиномы очень прилично изображают такие "дергунчики" очень малым числом членов и не вызывают побочных эффектов в виде ряби.

 

И почему посмотреть-то нельзя - я по временному ряду кита от стаи селёдки все равно не отличу, а любопытство заело...

Да там и смотреть нечего. Шум типа белого с какой-то шириной амплитуды. Да еще медленный тренд в течение суток. Амплитуда шума, кстати, тоже от времени суток зависима - днем она больше, ночью меньше. Ну а если где-то ... э ... типа зелетрясение произойдет али подземный ядерный взрыв :), но из того шума сразу вылезут импульсы, на остальной шум ничуточки не похожие. Т.е. сейсмические данные - очень близкий аналог моей задачи, хотя в моем случае "прослушивают" не землю, а живой организм.

 

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

 

С этой целью я хочу попросить вас растолковать мне явление, природа которого мне не ясна. А именно речь идет о дискретных представлениях полиномиального базиса (впрочем, у гармонического базиса тоже могут быть аналогичные проблемы). А сама проблема вот в чем. В аналитическом виде (в виде формулы) дела с ортогональностью обстоит лучше некуда - считают интегралы и те показывают ортогональность базиса. Но вот я загоняю тот базис в матрицу (скажем полином каждой степени занимает в ней отдельный столбец), то ортогональности как ни бывало! А точнее говоря, как бы я аккуратно не "оцифровывала" тот полином (хоть интегралом на 1/N-ой части области), всякий раз возникает значительная ошибка, как в отношении ортогональности между собой полиномов разных степеней, так и единичной нормы в отношении каждого из полинома.

 

Казалось бы причина этой бяки очевидна - своей дискретизацией я нарушаю что-то такое, что в аналитическом виде было гладким, а у меня стало ступенчатым. И в пользу такого заключения вроде бы свидетельствует тот факт, что с увеличением длины векторов (числа разбиений) ортогональность повышается. Казалось бы на коротких (меньше 10-ти элементов) и речи не может быть о точной ортогональности. Однако можно легко убедиться в том, что дискретизация ортогональности не помеха. Простой пример - берем квадратную матрицу полного ранга и просим библиотечную функцию найти полный набор ее собственных векторов. И вот абсолютно в любом случае (!) ортогональность тех векторов гарантирована с вычислительной точностью. Тут и точные нули между скалярными произведениями различных векторов будут на месте и единичная норма каждого из них. Причем, не приблизительная, а совершенно точная! И при этом на малых (коротких) векторах ничуть не хуже, чем на длинных.

 

Выходит, что дискретное представление не мешает выбирать ортогональные базисы. Тогда почему, когда дело касается полиномов Чебышева и Лежандра, то начинает катавасия с точностью? Это плохие полиномы? Или я их плохо оцифровала? Или они вообще для этих целей непригодны из-за того, что во времена Лежандра и Чебышева компьютеров не было и численными методами они не баловались? Или все-таки в дискретном виде могут существовать "гладкие" ортогональные многочлены?

 

Особняком тут стоят ряды/многочлены Фурье. Почему непонятно, но вроде бы у их дискретизация никаких проблем не вызывает. Хоть из 8-ми точек будет массив, но фурьевый базис как-то ухитряется быть ортогональным. Впрочем, я не исследовала случай, когда число точек не 2ⁿ, а другое. И потому не стану на их счет что-то категорически утверждать.

 

Что касается вида самих полиномов Чебышева и Лежандра, но они мне вполне нравятся, но отсутствие ортогональности (в векторном смысле) делает их для моих целей неприемлемыми. Где-то слышала, что бывают дискретные преобразования в базисы Чебышева и Лежандра (например, "Discrete Chebyshev Transform" и "Discrete Legendre Transform"), но ничего подходящего пока не нашла. Хотя в отношении Лежандра наклёвывается вариант "DLOP" (аббревиатура от "Discrete Legendre Orthogonal Polynomials"), но как следует разобраться с ним не успела.

[Xenia 42]

В отношении дискретного представления полиномов Лежандра (DLOP = "Discrete Legendre Orthogonal Polynomials") у меня наметился прогресс. Помогли статьи, найденные мною на IEEE:

1) On the computation of discrete Legendre polynomial coefficients (1992)

2) Recursive computation of discrete Legendre polynomial coefficients (1996)

 

DLOPы, какой бы короткой размерности ни были, сохраняют ортогональность идеально! Правда нормы у этого базиса по построению не единичные, но это я сама виновата - не осилила вычисление весовых коэффициентов. Но нормализовать их ничего не стоит обычным способом - вычислить норму в лоб и разделить на нее. Главное, что кросс-матрица у них изначально строго диагональна, как и положено.

 

Тем не менее имеются определенные странности, отличающие DLOP от обычных полиномов Лежандра. И этот вопрос заслуживает отдельного опуса :).

 

Во-первых, обычные полиномы Лежандра (как и полиномы Чебышева) определены на отрезке [-1, +1], тогда как DLOP - исключительно на положительных значениях натурального аргумента, т.е. определены дискретно на векторе/массиве ровно этой длины, для которой они строятся. Это очень удобно тем, что избавляет от необходимости постоянно пересчитывать номера элементов в X-координату на отрезок [-1, +1] и обратно.

 

Во-вторых, центр симметрии у обычных полиномов Лежандра (как и у полиномов Чебышева) находится точно в начале координат (в середине отрезка [-1,1]), тогда как у DLOP тоже в середине, но уже в середине вектора! Например, если вектор у нас длиной в 15 элементов, но центр симметрии будет приходиться на 7-ой элемент массива.

 

В-третьих, обычные полиномы Лежандра (как и полиномы Чебышева) всегда пересекаются в точке (1,1), т.е. в правом верхнем углу графика. А DLOP пересекаются всегда в 1-ом элементе вектора (в языке C у него нулевой номер). Из-за этого некоторые полиномы нечетной степени (которые раньше упирались в точку(-1,-1)) оказались перевернутыми вверх ногами.

 

Но в целом - красотища! Смотрим картинки.

 

 

- это обычные полиномы Лежандра, срисованные из интернета.

 

А полиномы DLOP, построенные по 15 точкам (array[15]), я построила сама и попытаюсь вложить в конец этого сообщения. Несмотря некоторую "корявость" представления полиномов высших степеней (4 и 5) из-за недостаточной длины вектора, ортогональность здесь превосходна!

 

Что же касается практических применений, то тут показательна статья

Adaptive ECG Data Compression Using Discrete Legendre Transform (2002)

в которой DLOP использованы для компрессии ЭКГ-сигнала. Как выглядит кардиограмма человека, наверняка, все знают. И по внешнему виду схожа со многими сигналами, с которыми часто приходится иметь дело: типа того, что долго не было ничего, и вдруг вигля выскочила :). А применимость для компрессии подобного вида сигналов свидетельствует о том, что базис DLOP показал себя на деле даже лучше, чем требуется для целей сглаживания шумов.

[_Ivana 0]

Xenia спасибо, очень интересны ваши подробные описания собственных размышлений, поисков и нахождения решений! Да ещё и с тематическики ссылками!

 

ЗЫ я после попытки устроить подобное (хотя нельзя считать её неудачной) подумал что на этом форуме не принято совместное коллективное исследование и решение задач. Вы возвращаете мне надежду :)

 

ЗЗЫ по делу сказать пока ничего не могу, хотя я предполагал что должны существовать красивые дискретно задаваемые ортогональные базисы с хорошими показателями нормы и ортогональности - но вы их сами уже нашли.

 

ЗЗЗЫ Идея Вашего метода понятна, за исключением деталей, которые на первый взгляд несущественны (типа бежим ли окном смещаясь на один отсчет входных данных или сразу кусками почти на всю длину окна - для оставления кусков перекрытия)