Матричное поэлементное умножение со сдвинутыми столбцами

Question

Скажем, у меня есть два массива, A и B.

Умножение по элементам определяется следующим образом:

Я хочу сделать элемент-мудрое умножение сверточно-подобным образом, т. е. перемещать каждый столбец на один шаг вправо, например, столбец 1 теперь будет столбцом 2, а столбец 3 теперь столбцом 1. Это должно привести к массиву (2 на 3 на 3)Матрица 2х3 для всех 3 возможностей)

Answer 1

Мы можем объединить A с одним собственным срезом, а затем получить эти скользящие окна.Чтобы получить эти окна, мы можем использовать np.lib.stride_tricks.as_strided на основе scikit-image's view_as_windows.Затем умножьте эти окна на B для окончательного вывода. Больше информации об использовании as_strided на основе view_as_windows.

Следовательно, у нас будет одно векторизованное решение, например, так:

In [70]: from skimage.util.shape import view_as_windows

In [71]: A1 = np.concatenate((A,A[:,:-1]),axis=1)

In [74]: view_as_windows(A1,A.shape)[0]*B
Out[74]: 
array([[[1, 0, 3],
        [0, 0, 6]],

       [[2, 0, 1],
        [0, 0, 4]],

       [[3, 0, 2],
        [0, 0, 5]]])

Мы также можем использовать multi-cores с numexpr module для последнего шага broadcasted-multiplication, который должен быть лучше на больших массивах.Следовательно, для примера, это было бы -

In [53]: import numexpr as ne

In [54]: w = view_as_windows(A1,A.shape)[0]

In [55]: ne.evaluate('w*B')
Out[55]: 
array([[[1, 0, 3],
        [0, 0, 6]],

       [[2, 0, 1],
        [0, 0, 4]],

       [[3, 0, 2],
        [0, 0, 5]]])

Синхронизация на больших массивах, сравнивая предложенные два метода -

In [56]: A = np.random.rand(500,500)
    ...: B = np.random.rand(500,500)

In [57]: A1 = np.concatenate((A,A[:,:-1]),axis=1)
    ...: w = view_as_windows(A1,A.shape)[0]

In [58]: %timeit w*B
    ...: %timeit ne.evaluate('w*B')
1 loop, best of 3: 422 ms per loop
1 loop, best of 3: 228 ms per loop

---

Выжимая лучшееметод на основе шага

Если вы действительно выжмете лучшее из подхода на основе шага, используйте исходный метод np.lib.stride_tricks.as_strided, чтобы избежать перегрузки в работеview_as_windows -

def vaw_with_as_strided(A,B):
    A1 = np.concatenate((A,A[:,:-1]),axis=1)
    s0,s1 = A1.strides
    S = (A.shape[1],)+A.shape
    w = np.lib.stride_tricks.as_strided(A1,shape=S,strides=(s1,s0,s1))
    return w*B

По сравнению с базирующимся на @Paul Panzer's array-assignment кроссовер выглядит как массивы 19x19 -

In [33]: n = 18
    ...: A = np.random.rand(n,n)
    ...: B = np.random.rand(n,n)

In [34]: %timeit vaw_with_as_strided(A,B)
    ...: %timeit pp(A,B)
10000 loops, best of 3: 22.4 µs per loop
10000 loops, best of 3: 21.4 µs per loop

In [35]: n = 19
    ...: A = np.random.rand(n,n)
    ...: B = np.random.rand(n,n)

In [36]: %timeit vaw_with_as_strided(A,B)
    ...: %timeit pp(A,B)
10000 loops, best of 3: 24.5 µs per loop
10000 loops, best of 3: 24.5 µs per loop

Итак,для чего-то меньшего, чем 19x19, array-assignment кажется лучше, а для более крупного, чем пошаговое, должен быть путь.

Answer 2

Просто заметка о view_as_windows / as_strided.Как бы ни были хороши эти функции, полезно знать, что они имеют довольно выраженные постоянные накладные расходы.Вот сравнение между решением @ Divakar view_as_windows (vaw) и моим подходом, основанным на копировании.

Как вы можете видеть, vawне очень быстро работает с малыми и средними операндами и начинает светить только при размере массива 30x30.

Код:

from simple_benchmark import BenchmarkBuilder, MultiArgument
import numpy as np
from skimage.util.shape import view_as_windows

B = BenchmarkBuilder()

@B.add_function()
def vaw(A,B):
    A1 = np.concatenate((A,A[:,:-1]),axis=1)
    w = view_as_windows(A1,A.shape)[0]
    return w*B

@B.add_function()
def pp(A,B):
    m,n = A.shape
    aux = np.empty((n,m,2*n),A.dtype)
    AA = np.concatenate([A,A],1)
    aux.reshape(-1)[:-n].reshape(n,-1)[...] = AA.reshape(-1)[:-1]
    return aux[...,:n]*B

@B.add_arguments('array size')
def argument_provider():
    for exp in range(4, 16):
        dim_size = int(1.4**exp)
        a = np.random.rand(dim_size,dim_size)
        b = np.random.rand(dim_size,dim_size)
        yield dim_size, MultiArgument([a,b])

r = B.run()
r.plot()

import pylab
pylab.savefig('vaw.png')

Answer 3

Запустите цикл for для количества столбцов и используйте np.roll() вокруг оси = 1, чтобы сдвинуть столбцы и выполнить умножение матриц.

см. Принятый ответ inэта ссылка.Надеюсь, это поможет.

Answer 4

Я могу фактически заполнить массив с двух сторон двумя колонками (чтобы получить массив 2x5) и запустить conv2 с 'b' в качестве ядра, я думаю, что он более эффективен