使用scipy solve_ivp时遭遇大规模数组内存分配异常问题

使用scipy.integrate.solve_ivp的BDF方法，利用雅可比矩阵稀疏性后，计算时间随网格点数量线性增长，但当数组元素量略超1000万时程序崩溃。理论上所需内存不足0.5GB，系统有32GB内存，不应出现此问题。

错误信息：

RuntimeError: SUPERLU_MALLOC fails for buf in intCalloc() at line 159 in file /private/tmp/scipy-20230221-70667-mabqq4/scipy-1.10.1/scipy/sparse/linalg/_dsolve/SuperLU/SRC/memory.c

复现代码：

import os
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
import scipy
from scipy.integrate import solve_ivp
import time
from time import perf_counter
from scipy.sparse import dia_array


def time_integration( Num ):

    t_grid = np.linspace( start=0., stop=1., num=Num )
    data = [np.zeros( Num, dtype=np.float64 )]
    jac_sparsity = dia_array( (data, [0]), shape=(Num, Num), dtype=np.float64 )

    def fun_rhs( t, y ):
        dfdt = np.ones( len(t_grid) )
        return dfdt

    start_time = perf_counter()
    output = scipy.integrate.solve_ivp( fun_rhs, [ min(t_grid), max(t_grid) ], np.zeros( len(t_grid) ), 
                                   method='BDF', 
                                   t_eval=None, 
                                   jac_sparsity=jac_sparsity, 
                                   dense_output=False, 
                                   events=None,
                                   vectorized=False,
                                   args=None, atol=1.e-2, rtol=1.e-2 )
    stop_time = perf_counter()
    return stop_time - start_time


start_log = 1.
stop_log = 7.2   ## !!!!! if 6. (or even 7.) it works fine !!!!!
timeSteps_array = np.logspace(start=start_log, stop=stop_log, num=int( stop_log-start_log ) + 1, dtype=int)
print(f'array of time steps (len = {len(timeSteps_array)}): {timeSteps_array}')
print('')


start_time_script = perf_counter()
runtime_array = np.zeros( len(timeSteps_array) )
runtime_array = [time_integration( it ) for it in timeSteps_array]
stop_time_script = perf_counter()


print(f'runtime for each Nt = {runtime_array} sec')
print(f'the script took {stop_time_script - start_time_script} seconds to run.')
plt.loglog(timeSteps_array, runtime_array, lw=2., color='blue', label='num test')
plt.loglog(timeSteps_array, time_integration( timeSteps_array[1] ) * timeSteps_array/timeSteps_array[1], lw=1.5, ls='--', color='red', label='linear')
plt.xlabel('Total number of points', fontsize=13)
plt.ylabel('Time [sec]', fontsize=13)
plt.legend(frameon=False, fontsize=12)
plt.savefig('TimeScaling.pdf',format='pdf',bbox_inches='tight', dpi=200)
plt.show()

原因分析 #

问题根源在Scipy依赖的SuperLU库上：

• SuperLU在处理超大维度稀疏矩阵时，内部内存分配可能存在32位整数溢出问题——即使系统是64位，部分旧版本SuperLU仍用32位整数计算内存块大小，当矩阵维度超过1e7时，计算出的内存大小会溢出，导致malloc调用失败。
• 即便雅可比是对角稀疏结构，SuperLU仍会分配一些辅助数组用于求解，这些数组的大小计算可能触发溢出。

解决方案 #

1. 升级Scipy版本 #

Scipy 1.11及以上版本对稀疏线性代数后端做了优化，修复了部分SuperLU的内存分配问题。直接升级到最新稳定版，大概率能解决该问题。

2. 更换ODE求解器方法 #

将method='BDF'替换为method='Radau'——Radau方法同样支持稀疏雅可比，且内部使用的线性代数后端对大规模问题兼容性更好，能绕过SuperLU的限制。

3. 优化问题结构 #

你的示例中每个ODE分量完全独立（y'=1），可以直接利用numpy向量化运算，避免构建超大稀疏矩阵：

def fun_rhs(t, y):
    return np.ones_like(y)

这种方式不需要传递jac_sparsity，求解器会自动处理，内存开销大幅降低。

4. 调整求解器参数 #

适当增大atol和rtol（当前为1e-2），宽松的收敛条件会减少迭代次数和辅助内存的使用，可能避开内存分配瓶颈。

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