Интерполяция/аппроксимация по обратному средневзвешенному расстоянию

Inverse distance weighting - это способ многомерной интерполяции на нерегулярной сетке. Существует много вариаций метода, отличающихся как концептуальными, так и техническими аспектами. В пакете ALGLIB реализован локальный вариант метода, генерирующий C1-непрерывные интерполянты и обладающий умеренной трудоемкостью (O(N·logN) - время построения модели, O(logN) - время интерполяции).

Содержание

    Простейший алгоритм
    Модифицированный метод Шепарда: интерполяция на неравномерной сетке
    Модифицированный метод Шепарда: аппроксимация на неравномерной сетке
    Настройка алгоритма: выбор Nw и Nq
    Настройка алгоритма: выбор узловых функций
    Manual entries

Простейший алгоритм

Рассмотрим простейший вариант IDW-алгоритма, перед тем, как перейти к более сложному алгоритму, входящему в состав ALGLIB. Пусть x_i - это набор точек в пространстве размерности D, а y_i - значения функции в этих точках. Тогда IDW-интерполянт имеет вид:

В такой формулировке алгоритм обладает как достоинствами, так и недостатками. Достоинствами являются:

Крайняя простота реализации
Отсутствие параметров, нуждающихся в настройке
Работоспособность в пространстве любой размерности
Работоспособность на любой равномерной/неравномерной сетке. Алгоритм "не замечает" мультиколлинеарности данных, способен успешно работать, даже если все точки сетки вложены в подпространство меньшей размерности. При правильной реализации алгоритм успешно работает даже на сетках с совпадающими узлами!

Однако недостатки алгоритма очень серьезны:

Низкая скорость работы на больших объемах данных: вычисление интерполянта имеет трудоемкость O(N).
Алгоритм придает слишком много веса удаленным узлам. Их суммарный вес может оказаться больше, чем вес узлов, расположенных рядом с точкой интерполяции. Это особенно выражено в пространствах высокой размерности
Глобальность интерполяции сама по себе является проблемой, т.к. интерполянт становится более чувствителен даже к далеким выбросам.
В узлах интерполяции функция f(x) плоская, т.е. имеет нулевую производную

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

Модифицированный метод Шепарда: интерполяция на неравномерной сетке

В пакете ALGLIB реализована более утонченная версия алгоритма, являющаяся незначительной модификацией модифицированного метода Шепарда, предложенного Robert J. Renka ('Multivariate Interpolation of Large Sets of Scattered Data', 1988). Модифицированный интерполянт имеет вид:

Модифицированный метод Шепарда отличается от оригинального алгоритма тем, что:

Для интерполяции используется подмножество точек K мощностью N_w - множество N_w ближайших соседей точки x.
Весовые функции W_i(x) имеют более сложный вид. Теперь они обращаются в ноль на границе сферы, радиус которой равен радиусу множества N_w ближайших соседей точки x.
Вместо константных значений y_i во взвешенной сумме теперь участвуют узловые функции Q_i(x). Эти функции могут быть квадратичными, линейными или константными (по выбору пользователя). Q_i(x) получается в результате взвешенной аппроксимации по МНК множества N_q ближайших соседей точки x_i с ограничением Q_i(x_i)=y_i. Весовые коэффициенты задачи МНК выбираются по формуле (2).
Для поиска ближайших соседей используются kd-деревья, позволяющие находить ближайшего соседа за время O(logN).

Основными достоинствами модифицированного алгоритма являются:

Всего два нуждающихся в настройке параметра: N_w и N_q (см. ниже). Оптимальные значения находятся перебором, но существует простой рецепт для получения хороших значений.
Работоспособность в пространстве любой размерности.
Работоспособность на любой равномерной/неравномерной сетке. Алгоритм "не замечает" мультиколлинеарности данных, способен успешно работать, даже если все точки сетки вложены в подпространство меньшей размерности. Алгоритм успешно работает даже на сетках с совпадающими узлами.
f(x) почти всюду C1-непрерывна (см. ниже)
Приемлемое быстродействие на больших наборах точек (достигнутое благодаря использованию быстрых алгоритмов поиска ближайших соседей). Построение модели имеет трудоемкость O(N·logN), интерполяция - O(logN).
Отсутствие "плоских пятен" в окрестностях узлов (при использовании квадратичных или линейных узловых функций).
Локальность интерполяционной схемы. Значение f(x) зависит только от ближайших к точке x узлов, что решает ряд проблем и существенно улучшает качество интерполяционной схемы.
Возможность использовать алгоритм для решения задач регрессии, т.е. для работы с зашумленными данными (см. ниже).

Можно отметить следующие недостатки алгоритма:

Алгоритм превосходит по скорости оригинальный IDW-алгоритм только на больших наборах точек (порядка нескольких сотен) и в пространствах умеренной размерности (2-5). На меньших наборах точек сказываются накладные расходы, связанные с навигацией по дереву поиска, хотя быстродействие алгоритма по-прежнему остается приемлемым. В пространствах более высокой размерности начинается снижение быстродействия поиска по kd-дереву.
Хотя алгоритм работоспособен на любой сетке, быстродействие снижается, если точки вложены в подпространство меньшей размерности. Это связано с тем, что в таких случаях для решения задач МНК используется основанный на SVD солвер вместо более быстрого основанного на QR солвера.
В некоторых редких случаях f(x) может иметь разрывы. Разрыв возникает, если более чем N_w ближайших соседей точки x находятся на абсолютно одинаковом расстоянии от нее. В этом случае алгоритм поиска N_w ближайших соседей не сможет сделать однозначный выбор, и, в зависимости от оругления в ходе операций с плавающей точкой, в окрестностях x ближайшие соседи будут выбираться случайным образом. Однако f(x) в любом случае остается ограниченной, а с ростом N_w вероятность такой ситуации снижается (в частности, нижняя граница для N_w установлена такой, чтобы подобная ошибка в принципе не могла возникнуть на прямоугольной сетке).

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

Модифицированный метод Шепарда: аппроксимация на неравномерной сетке

Алгоритм естественным образом модифицируется для работы с зашумленными данными: достаточно лишь незначительно модифицировать задачу МНК, решением которой является узловая функция Q_i(x). Во-первых, исчезает ограничение Q_i(x_i)=y_i. Во-вторых, весовые коэффициенты, использующиеся для вычисления Q_i(x), полагаются равными 1. Таким образом, функция Q_i(x) в некотором смысле усредняет соседние с xi точки.

Замечание #1
Сами формулы (1), (2) и (3) остаются прежними. Все изменения относятся лишь к тому, как вычисляются коэффициенты A, b и c в формуле (3).

Полученная таким образом, f(x) уже не проходит через предписанные точки, т.е. мы имеем задачу аппроксимации. Этот алгоритм имеет смысл применять если - и только если - мы работает с зашумленными данными. В незашумленном случае аппроксимационная схема не имеет преимуществ перед интерполяционной.

Настройка алгоритма: выбор Nw и Nq

Параметры N_w и N_q отвечают за количество узлов, использующихся на стадии интерполяции (параметр N_w) и на стадии построения узловых функций (параметр N_q):

параметр N_w контролирует локальность алгоритма. Чем больше его значение, тем больше узлов используется для интерполяции (ниже скорость), и тем больше значение интерполянта зависит от значений функции в удаленных узлах (ниже качество). С другой стороны, слишком малое значение N_w сделает интерполянт более резким, слишком сильно зависящим от близлежащих узлов (например, N_w=1 приведет к nearest neighbor interpolation). Хорошее значение N_w обычно ненамного больше чем max(1.5·N_q,2^D+1).
параметр N_q отвечает за другой аспект локальности: он отвечает за количество узлов, использующихся для определения рельефа узловых функций. Хорошая узловая функция должна проходить через (x_i,y_i) и воспроизводить рельеф функции в окрестностях x_i. Слишком большое значение - и узловая функция станет глобальной, адаптирующейся к значениям функции в далеких узлах, в которых эта узловая функция не особенно будет нужна. Хорошее N_q обычно в 1.5-2 раза больше, чем количество свободных параметров у узловой функции: 1+D - у линейной функции, (D+2)·(D+1)/2 - у квадратичной функции.

Ситуация несколько меняется при решении задач аппроксимации. Смысл параметра N_w и рекомендации по его выбору остаются прежними. Однако параметр N_q имеет другой смысл: теперь он отвечает не только за локальность алгоритма, но и за усреднение шумов во входных данных. Чем больше уровень шума, тем больше должно быть значение N_q, чтобы справиться с ним.

Настройка алгоритма: выбор узловых функций

При решении задач интерполяции пользователь может выбирать между четырьмя типами узловых функций:

константной узловой функцией. Это функция, использующаяся в оригинальном алгоритме, и оставленная исключительно для тестирования.
линейной узловой функцией, коэффициенты которой определяются по МНК. Эта функция обеспечивает лучшее качество, чем константная функция.
квадратичной узловой функцией, коэффициенты которой определяются по МНК. Эта функция обеспечивает самое лучшее качество, при условии, что у вас достаточно данных для вычисления коэффициентов квадратичной аппроксимации. Если это не так, лучше использовать линейную узловую функцию.
"быстрой" линейной узловой функцией, коэффициенты которой определяются путем интерполяции градиента. Это альтернативный вариант, зарезервированный для случаев, когда скорость построения модели важнее, чем её качество. Коэффициенты этой функции не минимизируют сумму квадратов невязок, но часто дают приемлемые результаты при интерполяции. Однако этот алгоритм несколько менее устойчив к выбросам и шумам во входных данных, чем предыдущие два алгоритма, поэтому его следует использовать, только если вы заранее можете убедиться в его работоспособности.

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