HPC集群Python矩阵求逆触发段错误、内存错误问题求助

• 从日志记录的矩阵大小15454226432字节（约14.4G）可推算，加载的是约44000阶的float64类型方阵。
• 初始段错误（Segmentation fault）的两个核心诱因：
  • 使用的NumPy 1.21.6默认会调用多线程BLAS库加速计算，如果没有手动限制线程数，库会自动尝试使用节点所有CPU核心创建线程，容易触发多线程栈空间溢出，触发非法内存访问报错。
  • 代码中先将稀疏矩阵转为稠密矩阵、再创建同尺寸单位阵、再计算i-df生成临时稠密矩阵，这三步已经占用约43G内存，np.linalg.inv调用底层LAPACK例程做LU分解、求逆时还需要额外分配临时内存，若内存峰值触碰到申请上限也可能触发段错误。
• 更换为scipy.sparse.linalg.inv后报内存错误的核心原因：绝大多数稀疏矩阵的逆是完全稠密的，没有稀疏性，该函数本质是逐列求解n个线性方程组得到完整逆矩阵，过程中产生的稠密中间结果内存开销远高于稠密矩阵直接求逆，即使申请500G内存也容易触发OOM。

按优先级从高到低排列：

• 优先方案：不要显式求逆
  99%的科学计算场景下，求矩阵逆的最终目的是求解线性系统或者计算逆矩阵和向量/矩阵的乘积，完全不需要显式生成完整逆矩阵，这种方式内存占用可以降低1-2个数量级，速度快10倍以上，数值稳定性也更好。
  参考实现：

  import scipy.sparse
  import scipy.sparse.linalg as spla
  
  # 保持稀疏格式加载矩阵，不要转稠密
  df = scipy.sparse.load_npz("test_data.npz").tocsc()
  n = df.shape[0]
  A = scipy.sparse.eye(n, format='csc') - df
  # 提前做LU分解，后续求解多个右侧向量可以复用分解结果
  A_factorized = spla.splu(A)
  
  # 如果需要计算 y = A^{-1} @ x，直接调用求解即可，不需要生成逆矩阵
  # y = A_factorized.solve(x)
  

• 必选配置：如果必须显式生成完整稠密逆矩阵
  1. 固定线性代数库线程数，避免多线程溢出：在Slurm提交脚本中加载Python模块前添加如下环境变量，线程数设置为你申请的CPU核心数即可，不要超过--ntasks指定的值：

     export OMP_NUM_THREADS=8
     export MKL_NUM_THREADS=8
     export OPENBLAS_NUM_THREADS=8
     export VECLIB_MAXIMUM_THREADS=8
     export NUMEXPR_NUM_THREADS=8
     

  2. 不要使用scipy.sparse.linalg.inv，直接将矩阵转为稠密格式后用SciPy的稠密线性代数接口求逆，开启参数复用内存降低峰值开销：

     import scipy.linalg as la
     import scipy.sparse
     
     df = scipy.sparse.load_npz("test_data.npz")
     n = df.shape[0]
     # 直接计算A = I - df的稠密格式，避免额外生成单位阵的内存拷贝
     A_dense = (scipy.sparse.eye(n) - df).toarray()
     # overwrite_a=True表示复用输入矩阵内存存储结果，check_finite=False跳过NaN/INF检查节省时间
     A_inv = la.inv(A_dense, overwrite_a=True, check_finite=False)
     

  3. 内存申请预留足够峰值：上述实现的峰值内存约为80-100G，Slurm脚本中设置--mem=120G即可满足需求，不需要申请500G内存。
• 可选优化：如果矩阵是对称正定、对角占优等有特殊结构的矩阵，可以替换为对应结构的专用分解接口（比如Cholesky分解），进一步降低内存占用、提升计算速度。

内容的提问来源于stack exchange，提问作者Stücke