Блог

Я использую numpy для выполнения умножения матриц, и я не могу понять, как использовать numpy для умножения 3D-матриц.

Допустим, у меня есть матрица 3×3, a, и я умножаю ее на вектор 3×1, b. Это даст вектор 3×1, c.

Это делается в numpy с:

 # (3, 3) * (3, 1) -> (3, 1)
c = np.matmul(a, b)
 

Хорошо, теперь я хочу выполнить аналогичную операцию с 3D-матрицей, которая по существу представляет собой 2500 матриц 3×3. Прямо сейчас я делаю что-то с эффектом:

 # (2500, 3, 3) * (2500, 3, 1) -> list of (3, 1) vectors with length 2500
C = [np.matmul(a, b) for a, b in zip(A, B)]
 

который возвращает список (3, 1) векторов.

Я бы предпочел НЕ зацикливаться, а вместо этого полностью использовать векторизацию numpy и матричные / тензорные произведения. Есть ли какая-то операция, которую я могу выполнить…

 # (2500, 3, 3) * (2500, 3, 1) -> (2500, 3, 1)
np.<function>(A, B, <args>)
 

Я видел материал об использовании np.tensordot, но я не знаю, как устанавливать оси.

 np.tensordot(A, B, axes=???)
 

Для имеющегося у вас массива размерности 3 (или тензора ранга 3) вы можете использовать np.einsum doc для более сложных матричных умножений. В вашем конкретном случае вы можете использовать следующее

 >>> import numpy as np
>>> x = np.random.randint(0, 3, (3, 3, 3))  # shape (3, 3, 3)
>>> y = np.random.randint(0, 3, (3, 3, 3))  # shape (3, 3, 3)
>>> np.einsum('ijk,ikl->ijl', x, y)  # still shape (3, 3, 3)
 

В частности, einsum выражение 'ijk,ikl->ijl' означает, что для каждой i th-матрицы выполните обычное умножение матрицы jk,kl->jl и поместите результат в i th-запись в результирующем тензоре (ndarray). Более общая форма этого процесса может быть

 np.einsum('...jk,...kl->...jl', x, y)
 

где вы можете иметь произвольное количество измерений перед каждым имеющимся у вас тензором (ndarray).

Смотрите Следующий полный пример:

 >>> import numpy as np
>>> x = np.random.randint(0, 3, (3, 3, 3))  # shape (3, 3, 3)
>>> x
array([[[0, 0, 1],
        [2, 2, 1],
        [2, 1, 1]],

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

       [[2, 2, 2],
        [1, 1, 2],
        [0, 2, 2]]])
>>> y = np.random.randint(0, 3, (3, 3, 3))  # shape (3, 3, 3)
>>> y
array([[[0, 0, 1],
        [2, 1, 0],
        [0, 0, 2]],

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

       [[0, 2, 1],
        [0, 1, 0],
        [0, 2, 1]]])
>>> np.einsum('ijk,ikl->ijl', x, y)
array([[[ 0,  0,  2],
        [ 4,  2,  4],
        [ 2,  1,  4]],

       [[ 6,  8,  2],
        [ 8,  6,  3],
        [10,  8,  4]],

       [[ 0, 10,  4],
        [ 0,  7,  3],
        [ 0,  6,  2]]])
>>> np.einsum('...ij,...jk->...ik', x, y)
array([[[ 0,  0,  2],
        [ 4,  2,  4],
        [ 2,  1,  4]],

       [[ 6,  8,  2],
        [ 8,  6,  3],
        [10,  8,  4]],

       [[ 0, 10,  4],
        [ 0,  7,  3],
        [ 0,  6,  2]]])
 

np.matmul(A,B) работает просто отлично. Какие ошибки вы получили?

 In [263]: A,B = np.arange(24).reshape(2,3,4), np.arange(8).reshape(2,4,1)
 

Решение einsum :

 In [264]: np.einsum('ijk,ikl->ijl',A,B)
Out[264]: 
array([[[ 14],
        [ 38],
        [ 62]],

       [[302],
        [390],
        [478]]])
In [265]: _.shape
Out[265]: (2, 3, 1)
 

Решение matmul :

 In [266]: A@B
Out[266]: 
array([[[ 14],
        [ 38],
        [ 62]],

       [[302],
        [390],
        [478]]])
 

Ваш цикл:

 In [267]: [np.matmul(a, b) for a, b in zip(A, B)]
Out[267]: 
[array([[14],
        [38],
        [62]]),
 array([[302],
        [390],
        [478]])]
 

matmul Документы:

 If either argument is N-D, N > 2, it is treated as a stack of
matrices residing in the last two indexes and broadcast accordingly.