быстрый способ инвертировать или поставить матрицу kxnxn

Question

Существует ли быстрый способ вычисления инверсии матрицы kxnxn с использованием numpy (инверсия рассчитывается для каждого k-среза)?Другими словами, есть ли способ векторизация следующий код:

>>>from numpy.linalg import inv
>>>a-random(4*2*2).reshape(4,2,2)
>>>b=a.copy()
>>>for k in range(len(a)):
>>>    b[k,:,:] = inv(a[k,:,:])

Answer 1

Сначала о получении обратного. Я посмотрел как np.linalg.tensorinv, так и np.linalg.tensorsolve.

Я думаю, к сожалению tensorinv не даст вам то, что вы хотите. Массив должен быть квадратным. Это исключает то, что вы хотите сделать, потому что их определение квадрата: np.prod(a[:i]) == np.prod(a[i:]), где i равно 0, 1 или 2 (одна из осей массива в целом); это может быть дано как третий аргумент ind из tensorinv. Это означает, что если у вас есть общий массив из NxN матриц длины M, вам нужно иметь, например, (для i = 1) NxN == NxM, что в общем случае неверно (на самом деле это верно в вашем примере, но все равно не дает правильного ответа).

Теперь, может быть, что-то возможно с tensorsolve. Это, однако, потребует некоторых тяжелых строительных работ с матричным массивом a, прежде чем он будет передан в качестве первого аргумента tensorsolve. Потому что мы хотели бы, чтобы b было решением "уравнения матрицы-матрицы" a*b = 1 (где 1 - это массив единичных матриц), а 1 будет иметь такую ​​же форму, как a и * 1020. *, мы не можем просто предоставить a, который вы определили как первый аргумент tensorsolve. Скорее, это должен быть массив с формой (M, N, N, M, N, N) или (M, N, N, N, M, N) или (M, N, N, N, N, M ). Это необходимо, потому что tensorsolve будет умножаться на b по этим трем последним осям, а также суммироваться по ним, так что результат (второй аргумент функции) снова будет иметь форму (M, N, N).

Затем, во-вторых, о точечных продуктах (ваше название предполагает, что это тоже часть вашего вопроса). Это очень выполнимо. Два варианта.

Во-первых: это сообщение в блоге Джеймса Хенсмана дает несколько хороших советов.

Второе: мне лично больше нравится использовать np.einsum для ясности. E.g.:

a=np.random.random((7,2,2))
b=np.random.random((7,2,2))
np.einsum('ijk,ikl->ijl', a,b)

Это умножит матрицы на все 7 "матриц" в массивах a и b. Кажется, что он примерно в 2 раза медленнее, чем метод array из поста выше, но он все равно примерно в 70 раз быстрее, чем использование цикла for, как в вашем примере. На самом деле, для больших массивов (например, 10000 матриц 5x5) метод einsum кажется немного быстрее (не знаю почему).

Надеюсь, это поможет.