Отказ от ответственности: я в основном пишу этот пост с учетом синтаксических соображений и общего поведения. Я не знаком с аспектами памяти и процессора описанных методов, и я нацелен на этот ответ тех, у кого достаточно небольшие наборы данных, так что качество интерполяции может быть основным аспектом, который следует учитывать. Я знаю, что при работе с очень большими наборами данных более эффективные методы (а именно griddata
и Rbf
) могут оказаться невозможными.
Я собираюсь сравнить три вида методов многомерной интерполяции ( interp2d
/ splines griddata
и Rbf
). Я подвергну их двум видам задач интерполяции и двум видам базовых функций (точки, из которых должны быть интерполированы). Конкретные примеры продемонстрируют двумерную интерполяцию, но жизнеспособные методы применимы в произвольных измерениях. Каждый метод обеспечивает различные виды интерполяции; во всех случаях я буду использовать кубическую интерполяцию (или что-то близкое к 1 ). Важно отметить, что всякий раз, когда вы используете интерполяцию, вы вносите смещение по сравнению с вашими необработанными данными, а конкретные используемые методы влияют на артефакты, которые вы в конечном итоге получите. Всегда помните об этом и выполняйте интерполяцию ответственно.
Две задачи интерполяции будут
- повышающая дискретизация (входные данные находятся в прямоугольной сетке, выходные данные - в более плотной сетке)
- интерполяция разрозненных данных на регулярную сетку
Две функции (по домену [x,y] in [-1,1]x[-1,1]
) будут
- плавная и дружелюбная функция
cos(pi*x)*sin(pi*y)
:; диапазон в[-1, 1]
- злая (и, в частности, прерывистая) функция:
x*y/(x^2+y^2)
со значением 0,5 около начала координат; диапазон в[-0.5, 0.5]
Вот как они выглядят:
Сначала я продемонстрирую, как эти три метода ведут себя в этих четырех тестах, а затем подробно расскажу о синтаксисе всех трех. Если вы знаете, чего ожидать от метода, возможно, вы не захотите тратить свое время на изучение его синтаксиса (глядя на вас, interp2d
).
Данные испытаний
Для наглядности вот код, с помощью которого я сгенерировал входные данные. Хотя в этом конкретном случае я, очевидно, знаю функцию, лежащую в основе данных, я буду использовать ее только для генерации входных данных для методов интерполяции. Я использую numpy для удобства (и в основном для генерации данных), но одного scipy тоже будет достаточно.
import numpy as np
import scipy.interpolate as interp
# auxiliary function for mesh generation
def gimme_mesh(n):
minval = -1
maxval = 1
# produce an asymmetric shape in order to catch issues with transpositions
return np.meshgrid(np.linspace(minval,maxval,n), np.linspace(minval,maxval,n+1))
# set up underlying test functions, vectorized
def fun_smooth(x, y):
return np.cos(np.pi*x)*np.sin(np.pi*y)
def fun_evil(x, y):
# watch out for singular origin; function has no unique limit there
return np.where(x**2+y**2>1e-10, x*y/(x**2+y**2), 0.5)
# sparse input mesh, 6x7 in shape
N_sparse = 6
x_sparse,y_sparse = gimme_mesh(N_sparse)
z_sparse_smooth = fun_smooth(x_sparse, y_sparse)
z_sparse_evil = fun_evil(x_sparse, y_sparse)
# scattered input points, 10^2 altogether (shape (100,))
N_scattered = 10
x_scattered,y_scattered = np.random.rand(2,N_scattered**2)*2 - 1
z_scattered_smooth = fun_smooth(x_scattered, y_scattered)
z_scattered_evil = fun_evil(x_scattered, y_scattered)
# dense output mesh, 20x21 in shape
N_dense = 20
x_dense,y_dense = gimme_mesh(N_dense)
Плавная функция и повышающая дискретизация
Начнем с самого простого. Вот как работает повышающая дискретизация от сетки формы [6,7]
к одной из [20,21]
функций гладкого теста:
Несмотря на то, что это простая задача, между выходами уже есть тонкие различия. На первый взгляд, все три выхода разумны. Следует отметить две особенности, основанные на наших предварительных знаниях о базовой функции: средний случай griddata
искажает данные больше всего. Обратите внимание на y==-1
границу графика (ближайшую к x
метке): функция должна быть строго равна нулю (так y==-1
как это узловая линия для гладкой функции), но это не так для griddata
. Также обратите внимание на x==-1
границу графиков (сзади, слева): основная функция имеет локальный максимум (подразумевающий нулевой градиент около границы) в [-1, -0.5]
, но griddata
выходные данные явно показывают ненулевой градиент в этой области. Эффект неуловимый, но, тем не менее, это предвзятость. (ВерностьRbf
даже лучше с выбором по умолчанию радиальных функций, дублированных multiquadratic
.)
Злая функция и повышающая дискретизация
Немного сложнее выполнить передискретизацию нашей злой функции:
Между этими тремя методами начинают проявляться явные различия. Глядя на графики поверхности, можно увидеть явные ложные экстремумы, появляющиеся на выходе interp2d
(обратите внимание на два выступа с правой стороны нанесенной поверхности). Хотя griddata
и Rbf
кажется, производят аналогичные результаты на первый взгляд, последний , кажется, производит более глубокий минимум вблизи , [0.4, -0.4]
что отсутствует основной функции.
Однако есть один важный аспект, в котором Rbf
он намного лучше: он уважает симметрию базовой функции (что, конечно, также стало возможным благодаря симметрии выборочной сетки). Выходной сигнал griddata
нарушает симметрию точек выборки, что уже слабо видно в гладком случае.
Гладкая функция и разбросанные данные
Чаще всего требуется выполнить интерполяцию по разрозненным данным. По этой причине я ожидаю, что эти тесты будут более важными. Как показано выше, точки выборки были выбраны псевдоднородно в интересующей области. В реалистичных сценариях у вас может быть дополнительный шум при каждом измерении, и вам следует подумать, есть ли смысл для начала интерполировать необработанные данные.
Выход для гладкой функции:
Сейчас уже происходит что-то вроде шоу ужасов. Я вырезал выходные данные от interp2d
до промежуточных [-1, 1]
исключительно для построения графика, чтобы сохранить хотя бы минимальный объем информации. Понятно, что хотя некоторая основная форма присутствует, есть огромные шумные области, где метод полностью не работает. Во втором случае griddata
форма довольно хорошо воспроизводится, но обратите внимание на белые области на границе контурного графика. Это связано с тем, что он griddata
работает только внутри выпуклой оболочки точек входных данных (другими словами, он не выполняет никакой экстраполяции ). Я сохранил значение NaN по умолчанию для выходных точек, лежащих за пределами выпуклой оболочки. 2 С учетом этих характеристик, Rbf
кажется, работает лучше всего.
Злая функция и разрозненные данные
И тот момент, которого мы все ждали:
Это не удивительно, что я interp2d
сдаюсь. Фактически, во время обращения к interp2d
вам следует ожидать, что некоторые друзья будут RuntimeWarning
жаловаться на невозможность построения сплайна. Что касается двух других методов, Rbf
кажется, что они дают лучший результат, даже вблизи границ области, в которой результат экстраполируется.
Итак, позвольте мне сказать несколько слов о трех методах в порядке убывания предпочтения (так, чтобы худший вариант был наименее вероятен для чтения кем-либо).
scipy.interpolate.Rbf
Rbf
Класса означает «радиальные базисные функции». Честно говоря, я никогда не рассматривал этот подход, пока не начал исследования для этого поста, но я почти уверен, что буду использовать его в будущем.
Как и в методах на основе сплайнов (см. Ниже), использование происходит в два этапа: первый создает вызываемый Rbf
экземпляр класса на основе входных данных, а затем вызывает этот объект для заданной выходной сетки для получения интерполированного результата. Пример из теста плавного повышения частоты дискретизации:
import scipy.interpolate as interp
zfun_smooth_rbf = interp.Rbf(x_sparse, y_sparse, z_sparse_smooth, function='cubic', smooth=0) # default smooth=0 for interpolation
z_dense_smooth_rbf = zfun_smooth_rbf(x_dense, y_dense) # not really a function, but a callable class instance
Обратите внимание, что и входные, и выходные точки в этом случае были 2-мерными массивами, а выходные данные z_dense_smooth_rbf
имеют ту же форму, что x_dense
и y_dense
без каких-либо усилий. Также обратите внимание, что Rbf
для интерполяции поддерживаются произвольные размеры.
Так, scipy.interpolate.Rbf
- производит корректный вывод даже для сумасшедших входных данных
- поддерживает интерполяцию в более высоких измерениях
- экстраполирует за пределы выпуклой оболочки входных точек (конечно, экстраполяция - это всегда азартная игра, и вам вообще не следует на нее полагаться)
- создает интерполятор в качестве первого шага, поэтому его оценка в различных выходных точках требует меньших дополнительных усилий
- могут иметь выходные точки произвольной формы (в отличие от прямоугольных сеток, см. ниже)
- склонны к сохранению симметрии входных данных
- поддерживает несколько видов радиальных функций для ключевого слова
function
: multiquadric
, inverse
, gaussian
, linear
, cubic
, quintic
, thin_plate
и определяемые пользователем произвольно
scipy.interpolate.griddata
Мой бывший фаворит, griddata
это обычная рабочая лошадка для интерполяции в произвольных измерениях. Он не выполняет экстраполяцию, кроме установки одного предустановленного значения для точек за пределами выпуклой оболочки узловых точек, но, поскольку экстраполяция - очень непостоянная и опасная вещь, это не обязательно является недостатком. Пример использования:
z_dense_smooth_griddata = interp.griddata(np.array([x_sparse.ravel(),y_sparse.ravel()]).T,
z_sparse_smooth.ravel(),
(x_dense,y_dense), method='cubic') # default method is linear
Обратите внимание на немного неуклюжий синтаксис. Входные точки должны быть заданы в массиве формы [N, D]
в D
размерах. Для этого мы сначала должны сгладить наши двумерные массивы координат (используя ravel
), затем объединить массивы и транспонировать результат. Есть несколько способов сделать это, но все они кажутся громоздкими. Входные z
данные также должны быть сглажены. Когда дело доходит до точек вывода, у нас немного больше свободы: по какой-то причине они также могут быть указаны как кортежи многомерных массивов. Обратите внимание, что help
of griddata
вводит в заблуждение, так как предполагает, что то же самое верно и для входных точек (по крайней мере, для версии 0.17.0):
griddata(points, values, xi, method='linear', fill_value=nan, rescale=False)
Interpolate unstructured D-dimensional data.
Parameters
----------
points : ndarray of floats, shape (n, D)
Data point coordinates. Can either be an array of
shape (n, D), or a tuple of `ndim` arrays.
values : ndarray of float or complex, shape (n,)
Data values.
xi : ndarray of float, shape (M, D)
Points at which to interpolate data.
В двух словах, scipy.interpolate.griddata
- производит корректный вывод даже для сумасшедших входных данных
- поддерживает интерполяцию в более высоких измерениях
- не выполняет экстраполяцию, для выхода может быть установлено одно значение вне выпуклой оболочки точек входа (см.
fill_value
)
- вычисляет интерполированные значения за один вызов, поэтому проверка нескольких наборов выходных точек начинается с нуля
- могут иметь выходные точки произвольной формы
- поддерживает интерполяцию ближайшего соседа и линейную интерполяцию в произвольных измерениях, кубическую в 1d и 2d. Метод ближайшего соседа и линейная интерполяция используются
NearestNDInterpolator
и LinearNDInterpolator
под капотом соответственно. 1d кубическая интерполяция использует сплайн, 2d кубическая интерполяция использует CloughTocher2DInterpolator
для построения непрерывно дифференцируемого кусочно-кубического интерполятора.
- может нарушить симметрию входных данных
scipy.interpolate.interp2d
/scipy.interpolate.bisplrep
Единственная причина, по которой я обсуждаю interp2d
и его родственников, заключается в том, что у него обманчивое название, и люди, вероятно, попытаются его использовать. Предупреждение о спойлере: не используйте его (начиная с версии 0.17.0). Он уже более особенный, чем предыдущие предметы, в том смысле, что он специально используется для двумерной интерполяции, но я подозреваю, что это наиболее распространенный случай для многомерной интерполяции.
Что касается синтаксиса, то interp2d
он похож на Rbf
то, что сначала необходимо создать экземпляр интерполяции, который можно вызвать для предоставления фактических интерполированных значений. Однако есть одна загвоздка: выходные точки должны быть расположены на прямоугольной сетке, поэтому входы, входящие в вызов интерполятора, должны быть 1d векторами, которые охватывают выходную сетку, как если бы из numpy.meshgrid
:
# reminder: x_sparse and y_sparse are of shape [6, 7] from numpy.meshgrid
zfun_smooth_interp2d = interp.interp2d(x_sparse, y_sparse, z_sparse_smooth, kind='cubic') # default kind is 'linear'
# reminder: x_dense and y_dense are of shape [20, 21] from numpy.meshgrid
xvec = x_dense[0,:] # 1d array of unique x values, 20 elements
yvec = y_dense[:,0] # 1d array of unique y values, 21 elements
z_dense_smooth_interp2d = zfun_smooth_interp2d(xvec,yvec) # output is [20, 21]-shaped array
Одна из наиболее распространенных ошибок при использовании interp2d
- это размещение ваших полных 2-мерных мешей в вызове интерполяции, что приводит к взрывному потреблению памяти и, надеюсь, к поспешному использованию MemoryError
.
Теперь самая большая проблема в interp2d
том, что это часто не работает. Чтобы понять это, нужно заглянуть под капот. Оказывается, interp2d
это оболочка для функций нижнего уровня bisplrep
+ bisplev
, которые, в свою очередь, являются оболочками для подпрограмм FITPACK (написанных на Фортране). Эквивалентный вызов в предыдущем примере будет
kind = 'cubic'
if kind=='linear':
kx=ky=1
elif kind=='cubic':
kx=ky=3
elif kind=='quintic':
kx=ky=5
# bisplrep constructs a spline representation, bisplev evaluates the spline at given points
bisp_smooth = interp.bisplrep(x_sparse.ravel(),y_sparse.ravel(),z_sparse_smooth.ravel(),kx=kx,ky=ky,s=0)
z_dense_smooth_bisplrep = interp.bisplev(xvec,yvec,bisp_smooth).T # note the transpose
Теперь вот что interp2d
: (в scipy версии 0.17.0) есть хороший комментарийinterpolate/interpolate.py
для interp2d
:
if not rectangular_grid:
# TODO: surfit is really not meant for interpolation!
self.tck = fitpack.bisplrep(x, y, z, kx=kx, ky=ky, s=0.0)
и действительно interpolate/fitpack.py
, bisplrep
там есть некоторая настройка и в конечном итоге
tx, ty, c, o = _fitpack._surfit(x, y, z, w, xb, xe, yb, ye, kx, ky,
task, s, eps, tx, ty, nxest, nyest,
wrk, lwrk1, lwrk2)
И это все. Базовые процедуры на interp2d
самом деле не предназначены для выполнения интерполяции. Их может быть достаточно для данных с достаточно хорошим поведением, но в реальных обстоятельствах вы, вероятно, захотите использовать что-то еще.
В заключение, interpolate.interp2d
- может привести к артефактам даже с хорошо закаленными данными
- специально для двумерных задач (хотя есть ограниченное количество
interpn
точек ввода, определенных в сетке)
- выполняет экстраполяцию
- создает интерполятор в качестве первого шага, поэтому его оценка в различных выходных точках требует меньших дополнительных усилий
- может производить вывод только по прямоугольной сетке, для разбросанного вывода вам придется вызывать интерполятор в цикле
- поддерживает линейную, кубическую и пятую интерполяцию
- может нарушить симметрию входных данных
1 Я абсолютно уверен , что cubic
и linear
вид базисных функций , Rbf
не точно соответствуют к другим интерполяторов одного и того же имени.
2 Эти NaN также являются причиной того, почему график поверхности кажется таким странным: у matplotlib исторически возникали трудности с построением сложных трехмерных объектов с правильной информацией о глубине. Значения NaN в данных сбивают с толку средство визуализации, поэтому части поверхности, которые должны быть сзади, отображаются так, чтобы быть впереди. Это проблема с визуализацией, а не с интерполяцией.