Как стать автором
Обновить

Комментарии 30

Это другая задача. Я решал не задачу интерполяции, а задачу аппроксимации. Интерполяция проходит через все точки, аппроксимация — нет. В том и в другом есть преимущества. Каждому инструменту — свое применение.

Иногда прохождение через все точки нежелательно. Данные могут содержать ошибки измерений. Невыгодно подгонять модель под случайные погрешности. Это приводит к проблемам, известным под названием "Overfitting".

За ссылку спасибо. Тоже интересный метод. Может пригодиться в подходящих обстоятельствах.

Справедливо. Я чесгря никогда не мог запомнить разницу :)

Классная штука. А для 3D подобное есть, не знаете?

Не знаю, к сожалению.

Для многомерной интерполяции я бы вам порекомендовал метод, основанный на радиально-базисных функциях, в английском языке RBF. При должном подборе параметров почти всегда можно обеспечить монотонность интерполянта между узлами. Я использовал этот метод для двумерной интерполяции с нерегулярной сеткой (гексагональная, с отдельными «дырками»). Очень доволен результатами.

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

Большое спасибо, никогда ранее не слышал про RBF

Стандартная задача на поиск условного минимума. Обычно решается применением метода множителей Лагранжа...

Я тоже думал применить к данной задаче метод множителей Лагранжа, но есть подозрение, что это приведет к нелинейным уравнениям, которые потом будет трудно решать. Детально не проверял.

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

Спасибо, обнадёжили! Надо будет попробовать. Если уравнения и правда останутся линейными — то с помощью метода множителей Лагранжа можно будет эффективно решать и более общую задачу линейной аппроксимации произвольными функциями с условием прохождения через точки. Достойно отдельной статьи!
НЛО прилетело и опубликовало эту надпись здесь
Эта статья рассчитана на подготовленного читателя, который уже знает (в университете это изучается), что такое метод наименьших квадратов, и умеет его применять (обычно с помощью компьютера — есть соответствующие программные пакеты и библиотеки).

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

1) для условия прохождения через заданные точки (а также значения производных в некоторых из этих точек) просто составлялась система уравнений. После её решения оставались несколько свободных коэффициентов, которые и можно было было подбирать.

2) Метод наименьших квадратов — не наилучший, а на наилёгкий, потому как легко математически описывается через производную интеграла. Наилучший — это минимизация максимального отклонения, но математически описать такое отклонение довольно проблематично (в общем случае, в частном случае многочлены Чебышева могут помочь). Однако некоторые инструменты численной аппроксимации позволяют к нему приблизиться.

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

Наилучший — это минимизация максимального отклонения

Чтобы говорить категориям лучший/худший надо ввести соответствующие определения. Обычно критерием является соответствующая метрика в пространстве функций. Так вот от её выбора может быть лучшим как метод наименьших квадратов, так и что-то иное.

Благодарю за внимание к статье.

2) Метод наименьших квадратов — не наилучший, а на наилёгкий

Насчет «наилучшего» хорошо ответил Naf2000 выше. Наилёгкий — это верно. И это важное свойство, которым можно и нужно пользоваться. Чтобы ставить и решать задачи, которые трудно или невозможно решать при других критериях оптимальности, таких как упомянутый вами минимакс.

Я тоже люблю минимакс-аппроксимацию и использую при случае. Алгоритм Ремеза — наше всё!

аппроксимация многочленом — тоже не наилучший, а наилёгкий способ

Безусловно, это наилёгкий способ. Но иногда он является и наилучшим — если истинная зависимость, которую моделируют по экспериментальным данным, является многочленом. Кроме того, многие другие случаи приближения сводятся к многочленам.

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

Каждому инструменту — своё применение. Конечно, большие степени для МНК обычно не применяются (а если применяются — то это классический Overfitting). Проблема пульсаций больше характерна не для МНК с многочленами, а для интерполяции многочленом Лагранжа. Опять же, надо владеть арсеналом разных методов и применять в конкретной ситуации наиболее подходящий. «Серебряной пули», универсального метода на все случаи жизни, не существует. Я увлекался минимакс-фильтрами для обработки сигналов, пока жизнь не подкинула такие сигналы, для которых этот метод плохо подходит.
Большие степени для МНК обычно не применяются
Я применял — для аппроксимации нарисованной мышкой АЧХ. Но мой комментарий был скорее о том, что классический подход хорошо работает на небольшом наборе данных и степени не выше 5. А как только требуется увеличить точность и количество исходных данных — сюрприз. Конечно, решение МНК через матрицы частично решают эту проблему, но это уже совсем другой уровень сложности.
Я увлекался минимакс-фильтрами для обработки сигналов, пока жизнь не подкинула такие сигналы, для которых этот метод плохо подходит.
Какие, если не секрет? Проблему фильтров я решил поиском вида функции (в частотной области), которая даёт максимально быстрое затухание импульсной характеристики (примерно e-x2/x).
Я применял — для аппроксимации нарисованной мышкой АЧХ.

Можно подробнее? Для чего аппроксимировалась АЧХ? Синтез фильтров?
классический подход хорошо работает на небольшом наборе данных и степени не выше 5

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

Что значит «через матрицы»? Бывает разве без матриц?
Какие, если не секрет?

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

Также оказались важны следующие свойства, нехарактерные для высококачественных аудио-фильтров: 1) сохранение амплитуды (т.е. прохождение интерполяционной кривой через исходные точки); 2) отсутствие колебаний импульсной характеристики.
Проблему фильтров я решил поиском вида функции (в частотной области), которая даёт максимально быстрое затухание импульсной характеристики

Почему именно импульсной? Это же наверняка приводило к ухудшению приближения частотной характеристики. Какие у вас были фильтры, рекурсивные? Уравнения Юла-Уокера не пробовали использовать? В матлабе есть функция "yulewalk" — для нахождения коэффициентов рекурсивных фильтров по АЧХ.
Можно подробнее? Для чего аппроксимировалась АЧХ?
1) Есть класс эффектов под названием «свёрточный ревербератор», у которых эффект реверберации достигается свёрткой сигнала с импульсной характеристикой помещения. У них есть подкласс «эмулятор гитарного кабинета» с аналогичным принципом действия. Импульсные характеристики для них записываются на микрофон откликом на единичный импульс. Пользователь перебирает эти импульсы и выбирает с наиболее подходящим звучанием.

Но перебирать сотни и тысячи импульсных характеристик — сомнительное удовольствие особенно с учётом того, что импульсы для одного и того же кабинета сильно разнятся в зависимости от положения микрофона при записи. Было бы удобнее контролировать звучание более управляемым способом.

Меня также интересовала эмуляция звучания телевизора для эмулятора Zx Spectrum — а для телевизоров импульсных характеристик никто не снимал.

2) Другая задача — это сглаживание АЧХ, также полученной через микрофон, для коррекции динамиков/колонок/помещения.

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

Бывает разве без матриц?
У меня нет специального образования, с математикой в школе было всё плохо (а с матрицами — особенно), не было доступа к литературе, доступ к интернету был пределом мечтаний, а про МНК я узнал на кухне за кружкой чая от человека, который тоже имел о нём весьма поверхностное уравнение. Поэтому я просто решил решил в Mathcad-е уравнения и получил формулы для коэффициентов в явном виде, непосредственно через координаты точек. Для больших степеней использовал базис a*cos(k x)+b*sin(k x) итеративным образом, последовательно увеличивая k.

МНК также можно посчитать чисто аналитически, для непрерывной функции, а не на дискретных значениях.

Почему именно импульсной? Это же наверняка приводило к ухудшению приближения частотной характеристики. Какие у вас были фильтры, рекурсивные?
FIR-фильтры с линейной фазой и наоборот, чисто фазовые фильтры с линейной АЧХ (например, преобразование Гильберта или ЛЧМ-подобный сигнал).

Проблема в том, что при получении коэффициентов фильтра через обратное преобразование Фурье они довольно медленно затухают (классический пример — sinc для brickwall фильтра) и их ограничивают умножением на оконную функцию, что приводит к искажению АЧХ. В простых случаях можно просчитать что получится, в сложных — уже сложнее, фазовые фильтры с линейной АЧХ таким образом вообще сделать не получится. Я искал способ, чтобы после обратного преобразования Фурье импульсная характеристика естественном образом уходила в ноль (в рамках погрешности вычислений). Это избавляет и от мучительных размышлений, какую именно оконную функцию выбрать, и позволяет определить необходимое количество коэффициентов для каждого конкретного фильтра.
Спасибо за ответы. Значит, вы со звуком много работали! Мне методы обработки звука тоже интересны.
эмуляция звучания телевизора для эмулятора Zx Spectrum

Ого, необычный проект! Вы участвуете в создании эмуляторов ZX Spectrum?
недостаточно просто удалить — нужно ещё чем-то заполнить освободившееся пространство

Линейное предсказание не пробовали? Я пробовал, результаты понравились.
про МНК я узнал на кухне за кружкой чая

«Классический» МНК сводится к решению системы линейных уравнений, а это и есть матричная алгебра в чистом виде. Можно, конечно, решать и без привлечения понятия матриц, но все равно они там незримо присутствуют.
МНК также можно посчитать… для непрерывной функции

Да, можно. Для синтеза FIR-фильтров есть такой метод.

Похоже, по синтезу фильтров у нас с вами общие интересы, но, наверно, общение на тему фильтров лучше будет перенести в личку или на более подходящий форум, а то тут несколько оффтопик. Пишите!
Вы участвуете в создании эмуляторов ZX Spectrum?
Нет, я просто брал готовые проекты с открытым исходным кодом и внедрял туда свои фильтры для проверки концепции.
пример
К сожалению, они все довольно ресурсоёмкие и для полноценного использования нужна оптимизация и куча дополнительного времени.

Линейное предсказание не пробовали?
Идея было в том, чтобы в случае «разорванной» синусоиды произвольной частоты восстанавливать потерянный фрагмент как можно ближе к оригиналу.
типа такого
Не уверен, что линейное предсказаниие так может.
Проблема пульсаций больше характерна не для МНК с многочленами, а для интерполяции многочленом Лагранжа
Контр-пример с многочленом 10-ой степени:

И пульсации хороши видно, и увеличение степени многочлена проблемы не решит.

А вот аппроксимация рациональным многочленом, той же степени и с тем же количеством коэффициентов:


И по поводу выбора подходящего инструмента: я не видел ни одной статьи на хабре о рациональной аппроксимации, зато статей о МНК миллиард. Так что нет никакого выбора.

UPD: вариант с шумом

UPD: соврал, степень рационального многочлена здесь в 2 раза меньше — числителя 5, знаменателя 4.

Как-то всё сложно...

Метод наименьших квадратов же не прибит гвоздями к базису 1, x², x³,.... Это раз.

И два – можно построить через фиксированные точки кривую (например, многочлен минимальной степени). Получается, что нам останется аппроксимировать отклонение точек от этой кривой (со знаком). Т.е. функцию с несколькими фиксированными корнями. Берём базис из функций с такими же корнями (надо подсказывать, как он получается из многочленов?) и прогоняем стандартный метод наименьших квадратов. Задача решена.

Как вы перешли от полинома, проходящего через фиксированные точки (x[i]!=0, y[i]!=0), к полиному с фиксированными корнями?

Вычитанием, конечно.

Допустим, у нас есть три фиксированные точки (X0,Y0), (X1,Y1), (X2,Y2) и десять точек (x_i,y_i) для приближения. По трём фиксированным точкам строим полином (параболу в данном случае) P0(x) = ax^2+bx+c (интерполяционный полином Лагранжа).
Дальше для остальных точек (x_i, y_i) считаем y'_i = y_i - P0(x_i) и приближаем их МНК по базису
1*(x-X0)*(x-X1)*(x-X2)
x*(x-X0)*(x-X1)*(x-X2)
x^2*(x-X0)*(x-X1)*(x-X2)
...
(столько, сколько надо)

Получили полином P1(x) такой, что P1(X0)=P1(X1)=P1(X2)=0 и приближающий наши разности.

Наш ответ, очевидно, P(x) = P0(x) + P1(x).
P0(X_i)=Y_i, P1(X_i)=0 - с фиксированными точками всё хорошо
y_i -P(x_i) = y'_i + P0(x_i) - P(x_i) = y'_i - P1(x_i) - как раз то, что мы старались улучшить МНК, опять всё хорошо.

Спасибо, теперь понятно стало.

Именно этот метод, с несущественным отличием, и был реализован в статье. Раздел 2.4.
Вы фактически пересказали основное содержание статьи (разд. 2.4.). Красиво, однако! А мне надо ещё много учиться лаконичности.

Ну, я по анекдоту про математика пересказывал:

Дано: Дрова, спички, пустой чайник, река.
Задача: Вскипятить воду.
Физик - Взять чайник, пойти наполнить водой, разжечь костер, укрепить чайник, ждать.
Математик - то же самое.

Дано: горящие Дрова, спички, полный чайник, река.
Задача: Вскипятить воду.
Физик - укрепить чайник, ждать.
Математик - затушить огонь из чайника. Задача сведена к предыдущей.

Зарегистрируйтесь на Хабре, чтобы оставить комментарий

Публикации

Истории