Ответ 1

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

Я собираюсь сравнить три вида многомерных методов интерполяции (interp2d/splines, griddata и Rbf). Я буду подвергать их двум типам задач интерполяции и двум типам базовых функций (точки, из которых должны быть интерполированы). Конкретные примеры продемонстрируют двумерную интерполяцию, но жизнеспособные методы применимы в произвольных измерениях. Каждый метод обеспечивает различные виды интерполяции; во всех случаях я буду использовать кубическую интерполяцию (или что-то близкое ¹). Важно отметить, что всякий раз, когда вы используете интерполяцию, вы вводите предвзятость по сравнению с вашими необработанными данными, а используемые конкретные методы влияют на артефакты, которые вы получите. Всегда будьте в курсе этого и интерполируйте ответственно.

Две задачи интерполяции будут

• upsampling (входные данные находятся на прямоугольной сетке, выходные данные находятся на более плотной сетке)
• интерполяция рассеянных данных на регулярную сетку

Две функции (над областью [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.)

Функция зла и повышение дискретизации

Более сложная задача состоит в том, чтобы выполнить upsampling по нашей злой функции:

Четкие различия начинают проявляться среди трех методов. Глядя на поверхностные графики, на выходе из interp2d появляются явные паразитные экстремумы (обратите внимание на два горба с правой стороны от поверхности). В то время как griddata и Rbf, на первый взгляд, дают аналогичные результаты, последний, по-видимому, создает более глубокий минимум около [0.4, -0.4], который отсутствует в базовой функции.

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

Плавная функция и рассеянные данные

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

Выход для гладкой функции:

Теперь уже идет немного ужасное шоу. Я привязал вывод interp2d к [-1, 1] исключительно для построения, чтобы сохранить как минимум минимальное количество информации. Ясно, что в то время как некоторые из основных форм присутствуют, есть огромные шумные области, где метод полностью разрушается. Второй случай griddata воспроизводит форму довольно красиво, но обратите внимание на белые области на границе контура. Это связано с тем, что griddata работает только внутри выпуклой оболочки входных точек данных (другими словами, она не выполняет никакой экстраполяции). Я сохранил значение NaN по умолчанию для выходных точек, лежащих вне выпуклой оболочки. ² Учитывая эти особенности, 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

Обратите внимание, что в этом случае обе входные и выходные точки были 2d-массивами, а выход 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, является общей рабочей лошадкой для интерполяции в произвольных измерениях. Он не выполняет экстраполяцию, не устанавливая одно заданное значение для точек вне выпуклой оболочки узловых точек, но поскольку экстраполяция является очень изменчивой и опасной вещью, это не обязательно является con. Пример использования:

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. Для этого сначала нужно сгладить наши 2-координатные массивы (используя 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 и ее родственников, состоит в том, что она имеет обманчивое имя, и люди, вероятно, попытаются ее использовать. Предупреждение о спойлере: не используйте его (как в scipy версии 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 заключается в том, чтобы положить полные 2d-ячейки в интерполяционный вызов, что приводит к потреблению взрывной памяти и, надеюсь, к поспешному MemoryError.

Теперь самая большая проблема с interp2d заключается в том, что она часто не работает. Чтобы понять это, мы должны смотреть под капот. Оказывается, что interp2d является оберткой для функций нижнего уровня bisplrep + bisplev, которые в свою очередь являются обертками для подпрограмм FITPACK (написанных на языке Fortran). Эквивалентным вызовом предыдущего примера будет

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 для входных точек, определенных на сетке)
• выполняет экстраполяцию
• создает интерполятор в качестве первого шага, поэтому оценка его в разных точках вывода меньше усилий
• может производить только выход по прямоугольной сетке, для разбросанного вывода вам нужно будет вызвать интерполятор в цикле
• поддерживает линейную, кубическую и квинтовую интерполяцию
• может нарушить симметрию входных данных

¹ Я вполне уверен, что базовые функции linear и linear Rbf не соответствуют другим интерполяторам с тем же именем.
²Эти NaN также являются причиной того, почему поверхностный график кажется таким странным: matplotlib исторически испытывает трудности с построением сложных 3D-объектов с надлежащей информацией о глубине. Значения NaN в данных путают визуализатор, поэтому части поверхности, которые должны быть сзади, изображены спереди. Это проблема с визуализацией, а не с интерполяцией.