Julia中如何高效初始化常量填充的大型结构化稀疏矩阵

你之前使用的代码M = spdiagm((a-b), N, N) .+ b无法实现压缩存储的核心原因是：

• 常规稀疏矩阵（比如Julia默认的SparseMatrixCSC、SciPy的csr_matrix/csc_matrix）的隐式非存储项默认值固定为0
• 对稀疏矩阵执行全矩阵广播加b的操作时，所有原本隐式存储的0位置都会被判定为“非默认值”，触发全量元素显式存储，直接退化为稠密矩阵，10×10矩阵存储100个元素就是这个问题导致的。

这种对角元为常数、非对角元也为常数的矩阵属于结构化矩阵，完全不需要存储N²个元素，内存占用可以做到和矩阵阶数N无关的O(1)级别，根据使用场景选对应方案即可：

• 如果你只需要做矩阵向量乘、矩阵基本运算，不需要显式存储元素：直接用线性算子逻辑实现，不需要初始化矩阵。
  该矩阵和任意向量v的乘法可以直接拆分为b * sum(v) + (a - b) * v，运算复杂度为O(N)，比任何稀疏存储格式的运算速度都快，10000阶以上规模也没有内存压力。Julia生态可以直接写函数调用，Python生态可以用scipy.sparse.linalg.LinearOperator封装该逻辑，兼容所有稀疏矩阵运算接口。
• 如果你需要兼容常规矩阵的索引、运算接口，可以用惰性填充结构组合实现，不会触发稠密化：
  以Julia生态为例，用FillArrays.jl提供的惰性常量矩阵搭配单位矩阵即可：

  using FillArrays, LinearAlgebra
  # Fill(b, N, N) 只存储常量b和矩阵维度，不分配N²内存
  M = (a - b) * I + Fill(b, N, N)
  

  该结构取M[i,j]时会自动判断：i==j返回a，否则返回b，支持大部分线性代数运算，内存占用仅为几个标量，和N的大小无关。
• 如果你必须和标准稀疏矩阵格式做交互，不要对生成的对角稀疏矩阵做全量加b操作，所有涉及矩阵的运算手动拆分对角项、全b常数项分别计算即可，从根源上避免稠密化。

不要尝试修改标准稀疏矩阵的默认非存储值，这类修改会破坏所有底层稀疏运算的预设逻辑，极易引发计算错误。

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