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

【问题标题】：Creating a sparse matrix from lists of sub matrices (Python)从子矩阵列表创建稀疏矩阵（Python）
【发布时间】：2016-07-29 23:26:53
【问题描述】：

这是我的第一个 SO 问题。让我知道我是否可以问得更好:)

我正在尝试找到一种将稀疏矩阵列表拼接成更大块矩阵的方法。

我有 python 代码，可以逐个矩阵地生成方稀疏矩阵列表。在伪代码中：

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

由于每个单独的 Lx1、Lx2 等矩阵都是按顺序计算的，因此它们被附加到一个列表中——我找不到“即时”填充类似数组的对象的方法。

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

M += J[i,j] * [ Lxi *Lxj + Lyi*Lyj + Lzi*Lzj ]

对于 0

似乎通过（伪代码）一步计算所有笛卡尔积来对其进行矢量化：

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

会更快。但是，np.bmat、np.vstack、np.hstack 等选项似乎需要数组作为输入，而我有列表。

有没有办法有效地将三个矩阵列表拼接成一个块？或者，有没有办法一次生成一个稀疏矩阵数组，然后将它们 np.vstack 在一起？

参考：用于计算 n-spin NMR 模拟的哈密顿矩阵的类似 MATLAB 代码可在此处找到：

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

【问题讨论】：

欢迎来到 Stack Overflow！第一个问题问得好。 :-)
小矩阵从何而来，采用什么稀疏矩阵格式？
每个 Lxn 矩阵由这些基本矩阵的重复 kron 计算：code sigma_x = csc_matrix(np.matrix([[0, 1/2], [1/2, 0]]) ) sigma_y = csc_matrix(np.matrix([[0, -1j/2], [1j/2, 0]])) sigma_z = csc_matrix(np.matrix([[1/2, 0], [0, - 1/2]])) 单位 = csc_matrix(np.matrix([[1, 0], [0, 1]])) code
（抱歉回复中的格式：我是新人！:)）
每个 Lxn 矩阵都是由这些基本矩阵的重复克朗计算的：sigma_x = csc_matrix(np.matrix([[0, 1/2], [1/2, 0]]))、sigma_y = csc_matrix(np.matrix([[0, -1j/2], [1j/2, 0]]))、sigma_z = csc_matrix(np.matrix([[1/2, 0], [0, -1/2]]))、unit = csc_matrix(np.matrix([[1, 0], [0, 1]]))

标签： python numpy scipy linear-algebra sparse-matrix

【解决方案1】：

这是scipy.sparse.bmat：

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

【讨论】：

将所有内容融合到一个矩阵中。我实际上需要保持各个子矩阵 Lx[n]/Ly[n]/Lz[n] 之间的界限。例如，L[1,0] 将包含 [Lx1, Ly1, Lz1] x [Lx0, Ly0, Lz0]。然后，J[1,0] 中的数字将与存储在 L[1,0] 中的矩阵相乘。我想您可以将 J 值与巨型矩阵的切片相乘（从而重新创建子矩阵），但这似乎很笨拙。
@GeoffreySametz：你不能像你想的那样“保持边界”，尤其是稀疏矩阵，它们本质上是二维的。即使使用密集的数组，“保持边界”也会让事情变得非常笨拙。
谢谢。我已经开始实施切片方法，并且会看到它是如何进行的。
事实证明这是行不通的，因为求和矩阵的转置与块元素的转置不同：math.stackexchange.com/questions/246289/…
@GeoffreySametz：哦，确实。那不是你真正想要的L.T。我对您的伪代码的解释与您的不同。新版本的答案应该会产生正确的结果。

【解决方案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

【讨论】：

【解决方案3】：

通过不使用稀疏矩阵和使用数组数组，我能够将速度提高十倍。

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] 矩阵仍然是“捆绑”的，因此 Product 是一个 n x n 矩阵。这意味着 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

【讨论】：