如何在Python中高效创建N个稀疏矩阵并构建更大的稀疏块矩阵？

我需要创建一个MN×MN的稀疏块矩阵，每个子块均为M×M的稀疏矩阵：

• 主对角线上的N个块A(j)依赖于各自的块索引j
• 其余对角线上的块B、C为固定稀疏矩阵

我已用scipy.sparse.diags创建B、C，但不清楚如何高效规范地生成A(j)，也不知道如何用scipy.sparse.bmat处理任意N的情况来合并这些块。此外，我尝试用for循环存储N个稀疏对角矩阵时失败，示例代码如下：

A = np.zeros(5)
for j in range(0,5):
    A[j] = j*scipy.sparse.diag(np.ones(3),shape(3,3))

补充：我利用所有块均为对角矩阵的特性，通过scipy.sparse.spdiags直接构建了大矩阵，示例代码如下：

import numpy as np
import scipy.sparse

M = 3
N = 5
C = np.arange(0,M)
d0 = np.zeros(M*N)
d1 = np.zeros(M*N)
dm1 = np.zeros(M*N)

dmM = np.ones(M*N)
dM = np.ones(M*N)

def somefunction(j):
    return np.sin(j)* j**2

for j in range(0,N):
    d0[j*M:(j+1)*M] = somefunction(j) * C
    dm1[j*M:(j+1)*M] = j*C
    d1[j*M:(j+1)*M] = j* C

data = np.zeros([5,M*N])
data[0,:] = dmM
data[1,:] = dm1
data[2,:] = d0
data[3,:] = d1
data[4,:] = dM

Matrix = scipy.sparse.spdiags(data,[-M,-1,0,1,M], M*N, M*N)

1. 修复循环存储稀疏矩阵的错误 #

你之前用np.zeros(5)创建的是数值数组，无法存储稀疏矩阵对象，改用列表存储即可：

import numpy as np
import scipy.sparse

M = 3
N = 5
def somefunction(j):
    return np.sin(j) * j**2

A_list = []
for j in range(N):
    # 生成依赖j的M×M稀疏对角块A(j)
    diag_vals = somefunction(j) * np.arange(M)  # 对应原代码中的C数组
    A_j = scipy.sparse.diags(diag_vals, shape=(M, M), format='csr')
    A_list.append(A_j)

2. 用scipy.sparse.bmat构建任意N的块矩阵 #

实现步骤： #

• 先构建块矩阵的二维列表结构：
  • 主对角线位置填充A_list[j]
  • 次对角线（j行j-1列，j>0时）填充固定块B
  • 上对角线（j行j+1列，j<N-1时）填充固定块C
  • 其余位置填充None或空稀疏矩阵

示例代码：

# 定义固定的B、C块（可替换为任意M×M稀疏矩阵）
B = scipy.sparse.diags(np.ones(M), shape=(M, M), format='csr')
C = scipy.sparse.diags(np.ones(M)*2, shape=(M, M), format='csr')

# 初始化块矩阵结构
block_matrix = [[None for _ in range(N)] for _ in range(N)]

for j in range(N):
    # 主对角线块
    block_matrix[j][j] = A_list[j]
    # 次对角线块（j>0时）
    if j > 0:
        block_matrix[j][j-1] = B
    # 上对角线块（j<N-1时）
    if j < N-1:
        block_matrix[j][j+1] = C

# 合并为MN×MN的大稀疏矩阵
big_matrix = scipy.sparse.bmat(block_matrix, format='csr')

3. spdiags方法的优化 #

你利用块对角特性直接操作全局对角线的思路很高效，可优化代码可读性：

import numpy as np
import scipy.sparse

M = 3
N = 5
total_size = M * N

def somefunction(j):
    return np.sin(j) * j**2

# 生成各全局对角线的数据
# 主对角线（偏移量0）
d0 = np.concatenate([somefunction(j) * np.arange(M) for j in range(N)])
# 相邻元素对角线（偏移量1和-1）
d1 = np.concatenate([j * np.arange(M) for j in range(N)])
dm1 = d1.copy()
# 块间对角线（偏移量M和-M）
dM = np.ones(total_size)
dmM = np.ones(total_size)

# 整理数据与对应偏移量
data = [dmM, dm1, d0, d1, dM]
offsets = [-M, -1, 0, 1, M]

# 构建最终稀疏矩阵
big_matrix = scipy.sparse.spdiags(data, offsets, total_size, total_size, format='csr')

方法选择建议： #

• 若子块不全是对角矩阵，优先用bmat，支持任意稀疏子块的组合
• 若子块均为对角矩阵，用spdiags效率更高，内存占用更低

内容的提问来源于stack exchange，提问作者Aitorkiv