Мне нужно быстро вычислить матрицу, элементы которой получаются путем свертки фильтра с вектором для каждой строки, субдискретизации элементов результирующего вектора, а затем взятия скалярного произведения результата с другим вектором. В частности, я хочу вычислить
M = [conv(e_j, f)*P_i*v_i]_{i,j},
где i изменяется от 1 до n, а j изменяется от 1 до m. Здесь e_j — вектор-индикатор (строка) размера n с единицей только в столбце j, f — фильтр длины s, P_i — матрица (n+s-1) на k, которая выбирает соответствующие k элементов из свертка, а v_i — вектор-столбец длины k.
Требуется O (n * s) операций для вычисления каждой записи M, поэтому в целом O (n * s * n * m) для вычисления M. Для n = 6, m = 7, s = 3, одно ядро моего компьютера (8GLOPs) должен уметь вычислять M примерно за 0,094 микросекунды. Тем не менее, моя очень простая реализация cython, следуя примеру, приведенному в документации Cython, требует более 2 миллисекунд, чтобы вычислить пример с этими параметрами. Это разница примерно в 4 порядка!
Вот файл shar с реализацией Cython и тестовым кодом. Скопируйте и вставьте его в файл и запустите «bash ‹fname›» в чистом каталоге, чтобы получить код, затем запустите «bash ./test.sh», чтобы увидеть ужасную производительность.
cat > fastcalcM.pyx <<'EOF'
import numpy as np
cimport numpy as np
cimport cython
from scipy.signal import convolve
DTYPE=np.float32
ctypedef np.float32_t DTYPE_t
@cython.boundscheck(False)
def calcM(np.ndarray[DTYPE_t, ndim=1, negative_indices=False] filtertaps, int
n, int m, np.ndarray[np.int_t, ndim=2, negative_indices=False]
keep_indices, np.ndarray[DTYPE_t, ndim=2, negative_indices=False] V):
""" Computes a numrows-by-k matrix M whose entries satisfy
M_{i,k} = [conv(e_j, f)^T * P_i * v_i],
where v_i^T is the i-th row of V, and P_i samples the entries from
conv(e_j, f)^T indicated by the ith row of the keep_indices matrix """
cdef int k = keep_indices.shape[1]
cdef np.ndarray M = np.zeros((n, m), dtype=DTYPE)
cdef np.ndarray ej = np.zeros((m,), dtype=DTYPE)
cdef np.ndarray convolution
cdef int rowidx, colidx, kidx
for rowidx in range(n):
for colidx in range(m):
ej[colidx] = 1
convolution = convolve(ej, filtertaps, mode='full')
for kidx in range(k):
M[rowidx, colidx] += convolution[keep_indices[rowidx, kidx]] * V[rowidx, kidx]
ej[colidx] = 0
return M
EOF
#-----------------------------------------------------------------------------
cat > test_calcM.py << 'EOF'
import numpy as np
from fastcalcM import calcM
filtertaps = np.array([-1, 2, -1]).astype(np.float32)
n, m = 6, 7
keep_indices = np.array([[1, 3],
[4, 5],
[2, 2],
[5, 5],
[3, 4],
[4, 5]]).astype(np.int)
V = np.random.random_integers(-5, 5, size=(6, 2)).astype(np.float32)
print calcM(filtertaps, n, m, keep_indices, V)
EOF
#-----------------------------------------------------------------------------
cat > test.sh << 'EOF'
python setup.py build_ext --inplace
echo -e "%run test_calcM\n%timeit calcM(filtertaps, n, m, keep_indices, V)" > script.ipy
ipython script.ipy
EOF
#-----------------------------------------------------------------------------
cat > setup.py << 'EOF'
from distutils.core import setup
from Cython.Build import cythonize
import numpy
setup(
name="Fast convolutions",
include_dirs = [numpy.get_include()],
ext_modules = cythonize("fastcalcM.pyx")
)
EOF
Я подумал, что, возможно, причиной может быть вызов scipy convolve (я не уверен, что cython и scipy хорошо работают вместе), поэтому я реализовал свой собственный код свертки, аналогичный примеру в документации Cython, но это привело к тому, что общий код примерно в 10 раз медленнее.
Любые идеи о том, как приблизиться к теоретически возможной скорости, или причины, почему разница так велика?
