Множественная одномерная интерполяция в двумерном массиве без цикла
У меня есть массив np с ndim = 2. Мне нужно выполнить выборку интерполированных значений по одному измерению, и я хотел бы сделать это как можно эффективнее.
Я придумал это решение:
for i in range(my_array.shape[0]):
my_interp_array[i, :] = np.interp(sample_y , np.arange(array_size_y), my_array[i,:])
... где sample_y — некоторый неэквидистантный вектор выборки. Хоть это и дает желаемый результат, возможно, это очень неэффективное решение.
Я также попробовал scipy.interpolate.interp1d, как предлагали другие, вот так:
y = np.arange(array_size_y) # equidistant sampling vector
intf = interp1d(y, my_array) # interpolation function
my_interp_array = intf(np.tile(sample_y, (len(y), 1)))
... но на самом деле это занимает в 3000 раз больше времени, чем первая идея с циклом for, что меня смущает.
Есть ли у кого-нибудь идеи, как это улучшить и почему scipy-интерполяция для меня намного медленнее?
Заранее спасибо!
Обновлено: Я сделал воспроизводимый пример, как некоторые об этом просили. Следующий код дает мне результаты, которые показывают, что scipy.interpolate.interp1d работает в 800-1000 раз медленнее, чем метод с циклом for.
import time
import numpy as np
from scipy.interpolate import interp1d
size_x, size_y = 512, 512
# make dummy arrays
my_array = np.random.rand(size_x, size_y)
my_interp_array = np.zeros_like(my_array)
# make indices
iy = np.arange(size_y) # "normal" indices
sample_y = iy + (iy ** 2 / np.max(iy ** 2)) # indices from where to get the values
sample_y[-1] = sample_y[-2] # avoid index out of range error
# for loop method
start1 = time.time()
for i in range(size_y):
my_interp_array[i, :] = np.interp(sample_y, iy, my_array[i, :])
end1 = time.time()
# scipy interp1d method
start2 = time.time()
intf = interp1d(np.arange(size_y), my_array)
my_interp_array = intf(np.tile(sample_y, (size_y, 1)))
end2 = time.time()
# print times
print(f"interp1d is {(end2 - start2) / (end1 - start1)} times slower than the for loop")
Это менее плохо, чем коэффициент 3000, но, на мой взгляд, разница все равно заметна.
Вы можете выйти из цикла np.arange(array_size_y), поскольку он кажется постоянным и каждый раз пересчитывается.
Имейте в виду, что общие интерполяции, как правило, требуют больших затрат на вычисление, особенно если область массива велика. Если ваши массивы большие, может помочь использование Numba с несколькими потоками. Если домен большой, то вычисления в любом случае будут ограничены памятью (из-за задержки памяти или пропускной способности из-за накладных расходов на размер строки кэша).
В случае, когда interp1d запрашивают линейную интерполяцию двумерного массива, он использует следующий код: github.com/scipy/scipy/blob/v1.14.0/scipy/interpolate/… Этот код векторизован, поэтому в этом отношении он может быть быстрее, чем цикл Python, но он также может иметь худшую эффективность кэширования из-за нескольких больших промежуточных массивов. Однако я был бы удивлен, если бы это привело к разнице в 3000 раз. Подтверждаю воспроизводимый пример.
Спасибо за ваши ответы. @jared Я разместил пример в редактировании. Что вы думаете?
@Nick Odell Это может быть объяснением этого, поскольку массивы, которые я использую, довольно большие (~ 512x512).
Вам следует использовать scipy.interpolate.make_interp_spline следующим образом: f = make_interp_spline(iy, my_array, k=1, axis=1) (k=1 означает линейную интерполяцию, и вы хотите, чтобы axis=1 соответствовало интерполяции, которую вы выполняете в версии с циклом), а затем интерполировать как my_interp_array = f(sample_y). И будьте осторожны: в вашем текущем коде эти два my_interp_arrays не равны.
@jared Готово ;) . Кстати, interp1d выполняет линейную интерполяцию, а scipy.interpolate.make_interp_spline представляет собой интерполяцию на основе сплайна (в вычислительном отношении значительно более затратная, хотя накладные расходы могут быть не одинаковыми и не эквивалентными).
Спасибо за ваш ответ. scipy.interpolate.make_interp_spline на самом деле значительно быстрее, чем interp1d, но все же медленнее, чем цикл for. Похоже, что временная сложность методов scipy увеличивается экспоненциально с увеличением размеров массива, в то время как цикл for увеличивается линейно (по крайней мере, вдоль одной оси). Я думаю, что пока буду придерживаться цикла.
@JérômeRichard С k=1 они дают те же результаты. Я проверил это с помощью np.allclose.
Ответы 1
Вы можете использовать scipy.interpolate.make_interp_spline с k=1 (для линейной интерполяции) и axis=1 (чтобы соответствовать той же оси интерполяции, что и код цикла).
from scipy.interpolate import make_interp_spline
f = make_interp_spline(iy, my_array, k=1, axis=1)
my_interp_array = f(sample_y)
Сравнение времени:
def original(iy, my_array, sample_y):
my_interp_array = np.zeros_like(my_array)
for i in range(size_y):
my_interp_array[i, :] = np.interp(sample_y, iy, my_array[i, :])
return my_interp_array
def mine(iy, my_array, sample_y):
f = make_interp_spline(iy, my_array, k=1, axis=1)
my_interp_array = f(sample_y)
return my_interp_array
%timeit original(iy, my_array, sample_y)
3.38 ms ± 82.1 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
%timeit mine(iy, my_array, sample_y)
1.64 ms ± 20.1 µs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)
Таким образом, этот подход приводит к увеличению скорости чуть более чем в 2 раза.
Моя машина ускоряется в 2,3 раза. Интересно, что 40% времени проводится в make_interp_spline, так что, думаю, там это не оптимально. Каждая функция Numpy на моей машине занимает 9 мкс, что настолько мало, что накладные расходы Numpy становятся значительными.
Я получаю разные (неправильные?) результаты с iy = np.sort(np.random.rand(512)*512). Это нормально?
@JérômeRichard Хм, иногда это работает, а иногда нет (в зависимости от того, что iy получится). Я не уверен, почему это так.
Удалось воспроизвести ваши результаты. Спасибо!
AFAIK
scipy.interpolate.interp1d— это удобная функция, но неэффективная, и она работает аналогично циклу, который вы выполняете с точки зрения производительности (возможно, с дополнительными накладными расходами, поскольку AFAIK более гибкая — например, нелинейная интерполяция). На самом деле, в документации указано: «Этот класс считается устаревшим и больше не будет получать обновления. Это также может означать, что он будет удален в будущих версиях SciPy». Таким образом, первое решение на самом деле лучше.