压缩格式下块对角矩阵的高效线性代数

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, p)和(N, p, m)的两个矩阵中实现有效的矩阵乘积，得到形状为(N, n, m)的矩阵？
• 是否存在一种方法可以使用numpy内置的（可能更快）函数执行这些操作？像np.linalg.inv这样的函数让我想到numpy被设计支持带压缩格式的对角块矩阵。

- dPol

---

听起来你也可以使用稀疏表示的标准矩阵乘法来处理分块对角系统。 - Bort

你是指使用 dot 方法的含义是什么？ - dPol

是的，可以将张量(N,n,m)作为压缩格式进行处理，也可以将它们保留为二维矩阵(Nn,Nm)并使用稀疏表示的稀疏矩阵乘法或矩阵向量乘法。 - Bort

你有想到哪个稀疏矩阵类？ - dPol

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是一种非常多才多艺的构造，用于复杂的线性运算，并且通常在两个操作数上效率很高。其思想是以索引方式重写您的操作：您需要的是对于每个i,j,l，将M1[i,j,k]与M2[i,k,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 是一个“块向量”，则结果将是一个块向量。

- Andras Deak -- Слава Україні

Answer 2

自Python 3.5及以上版本，可以使用矩阵乘法运算符@（numpy.matmul）简化上述示例，该运算符将此情况视为驻留在最后两个索引中的一堆矩阵并相应地进行广播：

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 = M1 @ M2  # similar to np.matmul(M1, M2)

all([np.allclose(np.dot(M1[k,...],M2[k,...]),Mprod[k,...]) for k in range(N)])
#True

- Hash04