scipy.linalg.block_diag 与 scipy.sparse.block_diag 在效率方面

scipy.linalg.block_diag vs scipy.sparse.block_diag in terms of efficiency

我对 scipy 构建分块对角矩阵的方式有疑问。我原以为创建一个稀疏块对角矩阵会比创建一个密集矩阵更快更有效(因为稀疏压缩)。但事实并非如此(或者我可能使用了一些低效的方法):

from timeit import default_timer as timer

import numpy as np
from scipy.sparse import block_diag as bd_sp
from scipy.linalg import block_diag as bd_la

m = [np.identity(1)] * 10000

before = timer()
res = bd_sp(m)
timer()-before
#takes 33.79 secs

before = timer()
res = bd_la(*m)
timer()-before
#takes 0.069 secs

我错过了什么?预先感谢您的回复。

In [625]: [np.identity(1)*i for i in range(1,5)]
Out[625]: [array([[1.]]), array([[2.]]), array([[3.]]), array([[4.]])]
In [626]: sparse.block_diag(_)
Out[626]: 
<4x4 sparse matrix of type '<class 'numpy.float64'>'
    with 4 stored elements in COOrdinate format>
In [627]: _.A
Out[627]: 
array([[1., 0., 0., 0.],
       [0., 2., 0., 0.],
       [0., 0., 3., 0.],
       [0., 0., 0., 4.]])

block_diag 使用 bmat 连接元素。 bmat 从所有元素中生成 coo 矩阵,并将它们的属性与偏移量组合,并生成一个新的 coo 矩阵。代码可读Python.

构建您自己的 data, row, col 数组可能更有效。 block_diag 很方便,适合合并几个大矩阵,但合并很多小矩阵时效率不高。

linalg 函数也是 Python(而且很短)。如果创建一个正确形状的 out 数组,并插入带有切片索引的块。这是一种高效的密集阵列解决方案。大部分艰苦的工作都是在编译后的 numpy 代码中完成的。

在进行矩阵乘法(和相关的 linalg 求解器)时,稀疏矩阵可以更快。对于大多数其他操作,包括初始化,它们比等效的密集代码慢。当问题太大时,它们也很有价值。