从子矩阵列表创建稀疏矩阵（Python）

Question

这是我的第一个问题。如果我能问得更好，请告诉我：）

我试图找到一种方法，把稀疏矩阵的列表拼接成一个更大的块矩阵。在

我有python代码，可以生成一个接一个矩阵的平方稀疏矩阵列表。在伪代码中：

Lx = [Lx1, Lx1, ... Lxn]
Ly = [Ly1, Ly2, ... Lyn]
Lz = [Lz1, Lz2, ... Lzn]   

由于每个单独的Lx1、Lx2等矩阵都是按顺序计算的，所以它们被附加到一个列表中——我无法找到一种“动态”填充类似数组的对象的方法。在

我正在优化速度，瓶颈是逐项计算笛卡尔积，类似于伪代码：

^{pr2}$

对于0<；=i的所有组合，j<；=n。（j是n维数字方阵）。在

通过一步计算所有笛卡尔积（伪码），可以将其矢量化：

L = [ [Lx1, Lx2, ...Lxn],
      [Ly1, Ly2, ...Lyn],
      [Lz1, Lz2, ...Lzn] ]
product = L.T * L

会更快。但是，诸如np.bmat公司, np.vstack公司, np.hstack公司似乎需要数组作为输入，而我有列表。在

有没有一种将矩阵拼接成三个矩阵的有效方法？或者，有没有办法一次生成一个稀疏矩阵数组，然后np.vstack公司他们在一起？在

参考文献：类似的MATLAB代码，用于计算n自旋NMR模拟的哈密顿矩阵，可以在这里找到：

http://spindynamics.org/Spin-Dynamics---Part-II---Lecture-06.php

Answer 1

这是^{}：

L = scipy.sparse.bmat([Lx, Ly, Lz], format='csc')
LT = scipy.sparse.bmat(zip(Lx, Ly, Lz), format='csr') # Not equivalent to L.T
product = LT * L

Answer 2

我有一个“矢量化”的解决方案，但它的速度几乎是原始代码的两倍。根据kernprof测试，上面显示的瓶颈和下面最后一行中显示的最终点积都占用大约95%的计算时间。在

    # Create the matrix of column vectors from these lists
L_column = bmat([Lx, Ly, Lz], format='csc')
# Create the matrix of row vectors (via a transpose of matrix with
# transposed blocks)
Lx_trans = [x.T for x in Lx]
Ly_trans = [y.T for y in Ly]
Lz_trans = [z.T for z in Lz]
L_row = bmat([Lx_trans, Ly_trans, Lz_trans], format='csr').T
product = L_row * L_column

Answer 3

使用稀疏矩阵和数组，我可以将速度提高10倍。在

Lx = np.empty((1, nspins), dtype='object') Ly = np.empty((1, nspins), dtype='object') Lz = np.empty((1, nspins), dtype='object')

在生成这些数组时，这些数组由单独的Lx数组（以前是稀疏矩阵）填充。使用数组结构可以使转置和笛卡尔积按需要执行：

Lcol = np.vstack((Lx, Ly, Lz)).real Lrow = Lcol.T # As opposed to sparse version of code, this works! Lproduct = np.dot(Lrow, Lcol)

单个Lx[n]矩阵仍然是“捆绑”的，因此乘积是nxn矩阵。这意味着n x n J数组与Lproduct进行就地乘法运算：

scalars = np.multiply(J, Lproduct)

然后将每个矩阵元素加到最终的哈密顿矩阵上：

for n in range(nspins): for m in range(nspins): M += scalars[n, k].real