Есть ли способ эффективно инвертировать массив матриц с помощью numpy?

Оказывается, вы сгораете на два уровня ниже в коде numpy.linalg. Если вы посмотрите на numpy.linalg.inv, вы увидите, что это просто вызов numpy.linalg.solve(A, inv(A.shape[0]). Это приводит к воссозданию матрицы идентичности на каждой итерации вашего для цикла. Поскольку все ваши массивы имеют одинаковый размер, это пустая трата времени. Пропуск этого шага путем предварительного выделения идентификационной матрицы сокращает время примерно на 20% (fast_inverse). Мое тестирование показывает, что предварительное выделение массива или выделение это из списка результатов не имеет большого значения.

Посмотрите на один уровень глубже, и вы найдете вызов подпрограммы lapack, но он завернут в несколько проверок работоспособности. Если вы удалите все это и просто вызовете lapack в своем цикле for (поскольку вы уже знаете размеры своей матрицы и, возможно, знаете, что она настоящая, а не сложная), все будет работать НАМНОГО быстрее (обратите внимание, что я сделал свой массив больше):

import numpy as np
A = np.random.rand(1000,3,3)
def slow_inverse(A): 
    Ainv = np.zeros_like(A)

    for i in range(A.shape[0]):
        Ainv[i] = np.linalg.inv(A[i])
    return Ainv

def fast_inverse(A):
    identity = np.identity(A.shape[2], dtype=A.dtype)
    Ainv = np.zeros_like(A)

    for i in range(A.shape[0]):
        Ainv[i] = np.linalg.solve(A[i], identity)
    return Ainv

def fast_inverse2(A):
    identity = np.identity(A.shape[2], dtype=A.dtype)

    return array([np.linalg.solve(x, identity) for x in A])

from numpy.linalg import lapack_lite
lapack_routine = lapack_lite.dgesv
# Looking one step deeper, we see that solve performs many sanity checks.  
# Stripping these, we have:
def faster_inverse(A):
    b = np.identity(A.shape[2], dtype=A.dtype)

    n_eq = A.shape[1]
    n_rhs = A.shape[2]
    pivots = zeros(n_eq, np.intc)
    identity  = np.eye(n_eq)
    def lapack_inverse(a):
        b = np.copy(identity)
        pivots = zeros(n_eq, np.intc)
        results = lapack_lite.dgesv(n_eq, n_rhs, a, n_eq, pivots, b, n_eq, 0)
        if results['info'] > 0:
            raise LinAlgError('Singular matrix')
        return b

    return array([lapack_inverse(a) for a in A])


%timeit -n 20 aI11 = slow_inverse(A)
%timeit -n 20 aI12 = fast_inverse(A)
%timeit -n 20 aI13 = fast_inverse2(A)
%timeit -n 20 aI14 = faster_inverse(A)

Результаты впечатляют:

20 loops, best of 3: 45.1 ms per loop
20 loops, best of 3: 38.1 ms per loop
20 loops, best of 3: 38.9 ms per loop
20 loops, best of 3: 13.8 ms per loop

РЕДАКТИРОВАТЬ: я недостаточно внимательно рассмотрел, что возвращается при решении. Получается, что матрица 'b' перезаписывается и в конце содержит результат. Этот код теперь дает согласованные результаты.

Несколько вещей изменилось с тех пор, как этот вопрос был задан и на него был дан ответ, и теперь numpy.linalg.inv поддерживает многомерные массивы, обрабатывая их как стопки матриц с последними индексами матриц (другими словами, массивы формы (...,M,N,N)). Похоже, это было появился в numpy 1.8.0. Неудивительно, что это лучший вариант с точки зрения производительности:

import numpy as np

A = np.random.rand(3,3,1000)

def slow_inverse(A):
    """Looping solution for comparison"""
    Ainv = np.zeros_like(A)

    for i in range(A.shape[-1]):
        Ainv[...,i] = np.linalg.inv(A[...,i])
    return Ainv

def direct_inverse(A):
    """Compute the inverse of matrices in an array of shape (N,N,M)"""
    return np.linalg.inv(A.transpose(2,0,1)).transpose(1,2,0)

Обратите внимание на два транспонирования в последней функции: ввод формы (N,N,M) должен быть транспонирован в форму (M,N,N) для работы np.linalg.inv, затем результат должен быть преобразован обратно в форму (M,N,N).

Результаты проверки и времени с использованием IPython на python 3.6 и numpy 1.14.0:

In [5]: np.allclose(slow_inverse(A),direct_inverse(A))
Out[5]: True

In [6]: %timeit slow_inverse(A)
19 ms ± 138 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

In [7]: %timeit direct_inverse(A)
1.3 ms ± 6.39 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

Циклы for действительно не обязательно намного медленнее, чем альтернативы, и в этом случае это вам не сильно поможет. Но вот предложение:

import numpy as np
A = np.random.rand(100,3,3) #this is to makes it 
                            #possible to index 
                            #the matrices as A[i]
Ainv = np.array(map(np.linalg.inv, A))

Сравнение времени этого решения с вашим решением дает небольшую, но заметную разницу:

# The for loop:
100 loops, best of 3: 6.38 ms per loop
# The map:
100 loops, best of 3: 5.81 ms per loop

Я попытался использовать подпрограмму numpy «vectorize» в надежде создать еще более чистое решение, но мне придется еще раз взглянуть на это. Изменение порядка в массиве A, вероятно, является наиболее значительным изменением, поскольку оно использует тот факт, что массивы numpy упорядочены по столбцам, и поэтому линейное считывание данных таким образом немного быстрее.

Вызовы Numpy-Blas не всегда являются самой быстрой возможностью

В задачах, где вам нужно вычислить множество обратных значений, собственных значений, скалярных произведений небольших матриц 3x3 или подобных случаев, numpy-MKL, который я использую, часто может быть значительно лучше.

Эти внешние подпрограммы Blas обычно создаются для задач с большими матрицами, для меньших вы можете написать стандартный алгоритм или посмотреть, например. Интел ИПП.

Также имейте в виду, что Numpy по умолчанию использует массивы C-упорядочения (последнее измерение изменяется быстрее всего). Для этого примера я взял код из Matrix inversion (3,3) python - жестко запрограммирован против numpy.linalg.inv и немного модифицировал его.

import numpy as np
import numba as nb
import time

@nb.njit(fastmath=True)
def inversion(m):
    minv=np.empty(m.shape,dtype=m.dtype)
    for i in range(m.shape[0]):
        determinant_inv = 1./(m[i,0]*m[i,4]*m[i,8] + m[i,3]*m[i,7]*m[i,2] + m[i,6]*m[i,1]*m[i,5] - m[i,0]*m[i,5]*m[i,7] - m[i,2]*m[i,4]*m[i,6] - m[i,1]*m[i,3]*m[i,8])
        minv[i,0]=(m[i,4]*m[i,8]-m[i,5]*m[i,7])*determinant_inv
        minv[i,1]=(m[i,2]*m[i,7]-m[i,1]*m[i,8])*determinant_inv
        minv[i,2]=(m[i,1]*m[i,5]-m[i,2]*m[i,4])*determinant_inv
        minv[i,3]=(m[i,5]*m[i,6]-m[i,3]*m[i,8])*determinant_inv
        minv[i,4]=(m[i,0]*m[i,8]-m[i,2]*m[i,6])*determinant_inv
        minv[i,5]=(m[i,2]*m[i,3]-m[i,0]*m[i,5])*determinant_inv
        minv[i,6]=(m[i,3]*m[i,7]-m[i,4]*m[i,6])*determinant_inv
        minv[i,7]=(m[i,1]*m[i,6]-m[i,0]*m[i,7])*determinant_inv
        minv[i,8]=(m[i,0]*m[i,4]-m[i,1]*m[i,3])*determinant_inv

    return minv

#I was to lazy to modify the code from the link above more thoroughly
def inversion_3x3(m):
    m_TMP=m.reshape(m.shape[0],9)
    minv=inversion(m_TMP)
    return minv.reshape(minv.shape[0],3,3)

#Testing
A = np.random.rand(1000000,3,3)

#Warmup to not measure compilation overhead on the first call
#You may also use @nb.njit(fastmath=True,cache=True) but this has also about 0.2s 
#overhead on fist call

Ainv = inversion_3x3(A)

t1=time.time()
Ainv = inversion_3x3(A)
print(time.time()-t1)

t1=time.time()
Ainv2 = np.linalg.inv(A)
print(time.time()-t1)

print(np.allclose(Ainv2,Ainv))

Производительность

np.linalg.inv: 0.36  s
inversion_3x3: 0.031 s