Комментарии 30
Существует интерполяция Steffen Spline - кривая проходит через все точки без дополнительных перегибов; даже есть реализация в GSL
Иногда прохождение через все точки нежелательно. Данные могут содержать ошибки измерений. Невыгодно подгонять модель под случайные погрешности. Это приводит к проблемам, известным под названием "Overfitting".
За ссылку спасибо. Тоже интересный метод. Может пригодиться в подходящих обстоятельствах.
Классная штука. А для 3D подобное есть, не знаете?
Не знаю, к сожалению.
Если же вас интересует многомерная аппроксимация — то метод наименьших квадратов обобщается на случай многомерных многочленов и, думаю, других базовых функций.
Стандартная задача на поиск условного минимума. Обычно решается применением метода множителей Лагранжа...
Вы к квадратичной форме добавите столько линейных условий (условие, что полином в данной точке равен фиксированному значению - линейное уравнение для коэффициентов), сколько точек у вас зафиксировано. Компоненты градиента модифицированной целевой функции - останутся линейными функциями от коэффициентов полинома и множителей Лагранжа.
В описаниях численных методов в литературе и научной периодике обычно приводится только теория и формулы. Опять же, подготовленному читателю именно это и нужно. Он сам потом сможет реализовать всё это в виде программы для компьютера.
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?Нет, я просто брал готовые проекты с открытым исходным кодом и внедрял туда свои фильтры для проверки концепции.
![](https://habrastorage.org/webt/ot/-r/4g/ot-r4g-_8coto5vk3dfdgcznxdu.jpeg)
Линейное предсказание не пробовали?Идея было в том, чтобы в случае «разорванной» синусоиды произвольной частоты восстанавливать потерянный фрагмент как можно ближе к оригиналу.
![](https://habrastorage.org/webt/2m/_z/h0/2m_zh0p0jj_m2ilyvjfdbs1rzww.png)
Проблема пульсаций больше характерна не для МНК с многочленами, а для интерполяции многочленом ЛагранжаКонтр-пример с многочленом 10-ой степени:
![](https://habrastorage.org/webt/yt/eb/we/ytebweb87htmex3xvegsdlqe8_s.png)
И пульсации хороши видно, и увеличение степени многочлена проблемы не решит.
А вот аппроксимация рациональным многочленом, той же степени и с тем же количеством коэффициентов:
![](https://habrastorage.org/webt/-u/gv/e8/-ugve8h4end6fhcewjdsiy90aye.png)
И по поводу выбора подходящего инструмента: я не видел ни одной статьи на хабре о рациональной аппроксимации, зато статей о МНК миллиард. Так что нет никакого выбора.
![](https://habrastorage.org/webt/jm/gr/l-/jmgrl-7wuwbpy7caxmagjswbbcm.png)
Как-то всё сложно...
Метод наименьших квадратов же не прибит гвоздями к базису 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) - как раз то, что мы старались улучшить МНК, опять всё хорошо.
Ну, я по анекдоту про математика пересказывал:
Дано: Дрова, спички, пустой чайник, река.
Задача: Вскипятить воду.
Физик - Взять чайник, пойти наполнить водой, разжечь костер, укрепить чайник, ждать.
Математик - то же самое.
Дано: горящие Дрова, спички, полный чайник, река.
Задача: Вскипятить воду.
Физик - укрепить чайник, ждать.
Математик - затушить огонь из чайника. Задача сведена к предыдущей.
Приближение многочленом с условием прохождения через точки