Эффективный способ смещения значений в одной матрице в соответствии со значениями в другой матрице

Question

Итак, я хочу сместить значения в matrix_a в соответствии со значениями в matrix_b. Поэтому, если значение в matrix_b в позиции 0,0 равно 1, то элемент в result_matrix в 0,0 должен быть элементом, который находится в 1,1 в matrix_a. У меня уже есть эта работа с использованием следующего кода:

import numpy as np

matrix_a = np.matrix([[1, 2,3], [4, 5, 6],[7, 8, 9]])
matrix_b = np.matrix([[1,1,0],[0,-1,0],[0,0,-1]])
result_matrix = np.zeros((3,3))


for x in range(matrix_b.shape[0]):
  for y in range(matrix_b.shape[1]):
    value = matrix_b.item(x,y)
    result_matrix[x][y]=matrix_a.item(x+value,y+value)

print(result_matrix)

, что приводит к:

[[5. 6. 3.]
 [4. 1. 6.]
 [7. 8. 5.]]

Сейчас это довольно медленно на больших матрицах, и у меня есть ощущение, что это может оптимизировать, используя одну из numpy или функции Сципи. Может кто-нибудь сказать мне, как это можно сделать более эффективно?

Answer 1

Использование np.indices

ix = np.indices(matrix_a.shape)
matrix_a[tuple(ix + np.array(matrix_b))]
Out[]: 
matrix([[5, 6, 3],
        [4, 1, 6],
        [7, 8, 5]])

В качестве совета, старайтесь избегать использования np.matrix - это действительно только для совместимости со старым MATLAB кодом и ломает много numpy функции. np.array работает так же хорошо в 99% случаев, а в остальное время np.matrix будет сбивать с толку основных numpy пользователей.

Answer 2

Вот один способ, когда integer-indexing генерируется из тех же итераторов, что и массивы открытого ранга, для получения индексов строк и столбцов для всех элементов -

I,J = np.ogrid[:matrix_b.shape[0],:matrix_b.shape[1]]
out = matrix_a[I+matrix_b, J+matrix_b]

Вывод для данного образца -

In [152]: out
Out[152]: 
matrix([[5, 6, 3],
        [4, 1, 6],
        [7, 8, 5]])

Синхронизация большого набора данных 5000x5000 -

In [142]: np.random.seed(0)
     ...: N = 5000 # matrix size
     ...: matrix_a = np.random.rand(N,N)
     ...: matrix_b = np.random.randint(0,N,matrix_a.shape)-matrix_a.shape[1]

# @Daniel F's soln
In [143]: %%timeit
     ...: ix = np.indices(matrix_a.shape)
     ...: matrix_a[tuple(ix + np.array(matrix_b))]
1.37 s ± 99.1 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

# Solution from this post
In [149]: %%timeit
     ...: I,J = np.ogrid[:matrix_b.shape[0],:matrix_b.shape[1]]
     ...: out = matrix_a[I+matrix_b, J+matrix_b]
1.17 s ± 3.21 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)