Julia中稀疏矩阵LU分解出现OutOfMemoryError的问题求助

需要对不同尺寸的稀疏矩阵执行LU分解作为GMRES预条件器，代码如下：

pc1 = ... # <15654×15654 SparseMatrixCSC{Float64, Int32} with 3327588 stored entries:
pc2 = ... # <31308×31308 SparseMatrixCSC{Float64, Int32} with 13310352 stored entries:

lu(pc1) # 运行正常
lu(pc2) # 抛出OutOfMemoryError

在190GB内存的超级计算机上重试仍报错，但将pc2转换为稠密矩阵后可成功执行LU分解，说明并非系统内存不足。堆栈跟踪显示错误来自UMFPACK库：

ERROR: OutOfMemoryError()
Stacktrace:
 [1] umferror(status::Int32)
   @ SparseArrays.UMFPACK /usr/local/Cellar/julia/1.9.0/share/julia/stdlib/v1.9/SparseArrays/src/solvers/umfpack.jl:111
 [2] macro expansion
   @ /usr/local/Cellar/julia/1.9.0/share/julia/stdlib/v1.9/SparseArrays/src/solvers/umfpack.jl:628 [inlined]
 [3] macro expansion
   @ ./lock.jl:267 [inlined]
 [4] umfpack_numeric!(U::SparseArrays.UMFPACK.UmfpackLU{Float64, Int32}; reuse_numeric::Bool, q::Nothing)
   @ SparseArrays.UMFPACK /usr/local/Cellar/julia/1.9.0/share/julia/stdlib/v1.9/SparseArrays/src/solvers/umfpack.jl:619
 [5] umfpack_numeric!
   @ /usr/local/Cellar/julia/1.9.0/share/julia/stdlib/v1.9/SparseArrays/src/solvers/umfpack.jl:618 [inlined]
 [6] #lu#7
   @ /usr/local/Cellar/julia/1.9.0/share/julia/stdlib/v1.9/SparseArrays/src/solvers/umfpack.jl:382 [inlined]
 [7] lu(S::SparseMatrixCSC{Float64, Int32})
   @ SparseArrays.UMFPACK /usr/local/Cellar/julia/1.9.0/share/julia/stdlib/v1.9/SparseArrays/src/solvers/umfpack.jl:378
 [8] top-level scope
   @ ./timing.jl:273 [inlined]
 [9] top-level scope
   @ ./REPL[53]:0

当前使用64位Julia版本。

错误根源在于默认使用的UMFPACK稀疏LU分解库：

• 对于Int32索引的矩阵，UMFPACK存在单个数组内存上限（约2GB），无法适配大规模矩阵的内存需求；
• UMFPACK的内存估算逻辑对某些高填充率的稀疏矩阵可能偏保守，导致提前触发内存不足错误；
• 稠密LU能成功执行，说明系统内存充足，问题出在UMFPACK的内存管理而非系统资源。

1. 切换为KLU稀疏LU求解器 #

KLU是SuiteSparse中针对大规模稀疏矩阵优化的求解器，内存管理更灵活，适合高维度稀疏矩阵：

using SuiteSparse.KLU
# 使用KLU算法执行LU分解
lu_result = lu(pc2, alg=KLU.KLUAlgorithm())

2. 改用Int64索引类型 #

将矩阵转换为Int64索引，突破UMFPACK在Int32下的内存寻址限制：

# 转换为Int64索引的稀疏矩阵
pc2_int64 = SparseMatrixCSC{Float64, Int64}(pc2)
# 执行LU分解
lu(pc2_int64)

3. 调整UMFPACK内存分配参数 #

手动增大UMFPACK的内存冗余因子，避免因估算不足导致的内存错误：

using SparseArrays.UMFPACK
# 初始化UMFPACK LU结构
umf_lu = UmfpackLU(pc2)
# 设置内存分配因子为2.0（默认值通常为1.2），增加冗余空间
umfset(umf_lu, :control, :alloc_factors, 2.0)
# 执行数值分解
umfpack_numeric!(umf_lu)

4. 分块/分布式LU分解（复杂场景） #

如果矩阵具有天然分块结构，可拆分后分别执行LU再组合；或使用分布式数组框架结合分布式求解器，适合超大规模矩阵场景。

优先尝试切换Int64索引和使用KLU后端，这两种操作简单且能覆盖绝大多数场景。若仍有问题，再调整UMFPACK的内存参数。

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