Numpy: эффективная линейная алгебра для блочно-диагональных матриц в сжатом формате

Question

У меня есть линейная система, в которой все матрицы являются диагональными. У них N блоки одинаковы по форме.

Матрицы хранятся в сжатом виде в виде массивов numpy с формой (N, n, m), а форма векторов - (N, m).

я в настоящее время реализована матрица-векторное произведение в

import numpy as np 

def mvdot(m, v): 
    return (m * np.expand_dims(v, -2)).sum(-1) 

Благодаря правилам вещания, если все блоки матрицы одинаковы я должен хранить его только один раз в массиве с формой (1, n, m): продукт с вектором (N, m) все еще дает правильный вектор (N, n).

Мои вопросы:

• как реализовать эффективную матрицу-матричный продукт, который дает матрицу с формой (N, n, m) из двух матриц с формами (N, n, p) и (N, p, m)?
• есть ли способ выполнить эти операции с помощью встроенной (возможно, более быстрой) функции numpy? Такие функции, как np.linalg.inv, заставляют меня думать, что numpy был разработан для поддержки этого сжатого формата для блочных диагональных матриц.

Answer 1

Если я правильно понимаю ваш вопрос, у вас есть два массива формы (N,n,p) и (N,p,m), соответственно, и их продукт должен иметь форму (N,n,m) где элемент [i,:,:] является матричным произведением M1[i,:,:] и M2[i,:,:]. Это может быть достигнуто с помощью numpy.einsum:

import numpy as np 
N = 7 
n,p,m = 3,4,5 
M1 = np.random.rand(N,n,p) 
M2 = np.random.rand(N,p,m) 
Mprod = np.einsum('ijk,ikl->ijl',M1,M2) 

# check if all the submatrices are what we expect 
all([np.allclose(np.dot(M1[k,...],M2[k,...]),Mprod[k,...]) for k in range(N)]) 
# True 

Numpy-х einsum невероятно универсален конструкция для сложных линейных операций, и это, как правило, довольно эффективным с двумя операндами. Идея состоит в том, чтобы переписать вашу операцию индексированным образом: вам нужно умножить M1[i,j,k] на M2[i,k,l] на каждые i,j,l и на сумму k. Это именно то, что делает вышеупомянутый вызов einsum: он сворачивает индекс k и выполняет необходимые продукты и задания по оставшимся размерам в данном порядке.

Матрицы-вектор продукт может быть сделан аналогичным образом:

M = np.random.rand(N,n,m) 
v = np.random.rand(N,m) 
Mvprod = np.einsum('ijk,ik->ij',M,v) 

Это возможно что numpy.dot может быть принужден с соответствующим транспонированием и уловками измерения непосредственно делать то, что вы хотите, но я не мог сделать та работа.

Оба указанных выше операций может быть сделано в том же вызове функции, позволяя неявное число измерений в пределах einsum:

def mvdot(M1,M2): 
    return np.einsum('ijk,ik...->ij...',M1,M2) 

Mprod = mvdot(M1,M2) 
Mvprod = mvdot(M,v) 

В случае, если входной аргумент M2 представляет собой блок-матрица, будет ведущим размер добавляется к результату, создавая матрицу блоков. В случае, если M2 является «блочным вектором», результатом будет вектор блока.