Numba编译场景下np.linalg.eigvalsh/eigvals异常触发不一致问题求助

我现在用numba编译一些计算量很大的代码，性能提升非常明显，这部分体验很棒！但最近我给计算逻辑加了个小改动——提取特征值，清掉缓存后开始测试，代码能正常编译运行，但时不时会毫无征兆地崩溃，抛出这个错误：

SystemError: <function _numba_unpickle at 0x0000027532E4EFC0> returned a result with an exception set

我的所有函数都是用下面的装饰器编译的：
适用普通场景用这个：

@njit(error_model="numpy", cache=True)

需要并行加速的场景用这个：

@njit(error_model="numpy", cache=True, parallel=True)

当前环境版本：python 3.12.9、numpy 2.1.3、numba 0.61.0，其他依赖是它们自动拉取的配套版本。

我已经把调试器设置为捕获所有抛出的异常，但帮助不大——看起来数据和代码本身都没问题：把编译函数里的逻辑拿到控制台跑完全没错误，甚至在调试器里调用触发过错误的那个编译函数，它也能正常运行。这个错误不是稳定复现的，完全随机。

调用触发过错误的编译函数时，会出现三种情况：

1. 返回正确结果
2. 返回完全无意义的垃圾值（比如xxxxx+303这样的超大数，大概率是特征值计算收敛失败了）
3. 抛出上面提到的异常

而且用相同数据重复测试时，错误出现的时机也随机：有时候调用1、2次就报错，有时候跑100次才会触发。偶尔还能抓到更具体的错误栈：

Traceback (most recent call last):
  File "C:\Users\jpmor\anaconda3\envs\new_base\Lib\site-packages\numba\np\linalg.py", line 841, in _check_finite_matrix
    raise np.linalg.LinAlgError(
numpy.linalg.LinAlgError: Array must not contain infs or NaNs.

The above exception was the direct cause of the following exception:

我的代码里只用到两次np.linalg相关操作：
一次是计算向量范数：

len = np.linalg.norm(vect[indx, :])

另一次就是特征值计算：

eigenvalues = np.linalg.eigvalsh(stress.reshape(3, 3))

（或者非对称场景下的np.linalg.eigvals）

从错误栈能看出来，_check_finite_matrix是np.linalg.eigvals的重载实现里调用的，所以问题肯定出在特征值计算上，但我完全搞不懂：相同的数据，为什么有时候报错有时候又能正常运行？

查Numpy文档知道，特征值计算收敛失败时会抛出LinAlgError，我怀疑这是收敛相关的问题，但奇怪的是：用纯Numpy的未编译代码跑相同数据，从来不会出现这个错误，只有用Numba编译后才会随机触发。

相关函数代码 #

核心应力计算函数

@njit(error_model="numpy", cache=True)
def compute_stress(num_elems, force, vect, weight_func, force_symmetry=False):
    results = np.zeros((1, 26))

    # 计算接触应力
    computed_stress = np.zeros((1, 9))
    for indx in range(num_elems):
        for a in range(3):
            for b in range(3):
                sid = 3 * a + b
                computed_stress[0, sid] += force[indx, a] * vect[indx, b] * weight_func[indx]

    if force_symmetry:
        # 对称化处理应力矩阵
        a = (computed_stress[0, 1] + computed_stress[0, 3]) / 2
        b = (computed_stress[0, 2] + computed_stress[0, 6]) / 2
        c = (computed_stress[0, 5] + computed_stress[0, 7]) / 2

        computed_stress[0, 1] = a
        computed_stress[0, 2] = b
        computed_stress[0, 5] = c
        computed_stress[0, 3] = a
        computed_stress[0, 6] = b
        computed_stress[0, 7] = c

        eigenvalues = np.linalg.eigvalsh(computed_stress.reshape(3, 3))

    else:
        eigenvalues = np.linalg.eigvals(computed_stress.reshape(3, 3))

    # 计算等效应力、偏应力等衍生量
    sigma_iso = np.trace(computed_stress.reshape(3, 3)) * np.eye(3) * 0.33333333333333333333
    sigma_dev = computed_stress.reshape(3, 3) - sigma_iso

    eigenvalues.sort()
    eigenvalues = eigenvalues[::-1]
    tau_max = (eigenvalues[0] - eigenvalues[-1]) / 2

    # 计算应力不变量
    J2 = 0.5 * np.trace(np.dot(sigma_dev, sigma_dev))

    # 其他常用应力指标
    mean_stress = np.mean(eigenvalues)  
    eqv_stress = np.sqrt(3.0 * J2)  

    stress_ratio = eigenvalues[0] / eigenvalues[-1]

    # 整理结果
    results[0, :9] = computed_stress
    results[0, 9:12] = eigenvalues
    results[0, 12] = tau_max
    results[0, 13] = mean_stress  
    results[0, 14] = eqv_stress  
    results[0, 15] = stress_ratio
    results[0, 16] = 0
    results[0, 17:] = sigma_dev.reshape(1, -1)

    return results

并行驱动函数

@njit(error_model="numpy", cache=True, parallel=True)
def compute_for_mesh(node_list, force_symmetry=False):

    node_results = np.zeros((node_list.shape[0], 39))

    for node_idx in prange(len(node_list)):
        
        # 此处省略其他前置计算逻辑
        # ...

        # 计算接触应力
        if num_elements > 0:
            
            # 此处省略调用compute_stress前的其他计算
            # ...

            stresses = compute_stress(num_elements, force, vect, weight_func, force_symmetry=force_symmetry)

            # 此处省略结果后续处理逻辑
            # ...

            node_results[node_idx, :26] = stresses

    return node_results

我的求助需求 #

我知道最简单的解决办法是不用Numba，但这样代码运行速度会慢到离谱，根本没法当天跑完。所以想请教：

1. 有没有和np.linalg.eigvals功能类似、且在Numba环境下更稳定的替代方案？
2. 如果自己写一个特征值计算函数，能不能做到和Numpy原生实现一样快？

备注：内容来源于stack exchange，提问作者jpmorr