从子矩阵列表创建稀疏矩阵 (Python)

Question

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

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

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

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

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

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

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

对于 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

这是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

Answer 2

我有一个 "vectorized" 解决方案，但它的速度几乎是原始代码的两倍。根据 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

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

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] 矩阵仍然是 "bundled"，因此 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

从子矩阵列表创建稀疏矩阵 (Python)

Creating a sparse matrix from lists of sub matrices (Python)

python

numpy

linear-algebra

scipy

sparse-matrix