Numpy大数运算问题：特征值与行列式计算异常排查

问题代码 #

import numpy as np
from numpy import exp


def foo(a,b):
    return np.array([[exp(-2*b*a), 0,-exp(b*a)/3 + exp(-2*b*a)/3],
                    [0, exp(-2*b*a), 0],
                    [-exp(b*a)/3 + exp(-2*b*a)/3,0, 10*exp(4*b*a)/9 - 2*exp(b*a)/9 + exp(-2*b*a)/9]])


a_s = np.linspace(0, 10, 100)
b_s = np.linspace(0.1, 5, 100)

store_det = np.zeros((len(a_s), len(b_s)))

for i, b in enumerate(b_s):
    for j, a in enumerate(a_s):
        store_det[i,j] = np.linalg.det(foo(a,b))
        
print(f'max: {np.max(store_det)} \t min: {np.min(store_det)}')


store_det_by_eigen = np.zeros((len(a_s), len(b_s)))

for i, b in enumerate(b_s):
    for j, a in enumerate(a_s):
        w,v = np.linalg.eig(foo(a,b))
        store_det_by_eigen[i,j] = w[0]*w[1]*w[2]

print(f'max: {np.max(store_det_by_eigen)} \t min: {np.min(store_det_by_eigen)}') 

问题描述 #

理论上该3x3矩阵对任意实数a、b，行列式值应为1，且特征值的乘积等于行列式。但当a=10、b=5时，特征值乘积不再为1，询问操作是否有误及解决办法。

原因分析 #

• 数值溢出与精度丢失：当a=10、b=5时，4*b*a=200，exp(200)的结果约为1e86，远超float64浮点数的有效范围，导致矩阵中对应元素变成inf，后续特征值计算的精度被严重破坏。
• 特征值相乘的误差放大：np.linalg.eig处理极端数值矩阵时，小特征值的计算精度会被大特征值掩盖，直接相乘后误差被进一步放大，偏离真实值1。

解决办法 #

1. 利用对数转换避免大数值直接相乘

利用行列式的对数等于特征值对数之和的性质，先计算特征值的对数和再指数还原，绕开超大数相乘的问题：

# 修改特征值乘积计算的循环部分
for i, b in enumerate(b_s):
    for j, a in enumerate(a_s):
        w, v = np.linalg.eig(foo(a,b))
        # 取绝对值避免极小负数值导致对数报错，再求和指数还原
        store_det_by_eigen[i,j] = np.exp(np.sum(np.log(np.abs(w))))

这种方法能有效规避大数值相乘的精度问题，计算结果会接近1，误差在浮点数允许范围内。

2. 矩阵归一化缩小数值差异

在生成矩阵时，将矩阵除以最大元素的量级，缩小元素间的数值差，减少精度丢失：

def foo(a,b):
    mat = np.array([[exp(-2*b*a), 0,-exp(b*a)/3 + exp(-2*b*a)/3],
                    [0, exp(-2*b*a), 0],
                    [-exp(b*a)/3 + exp(-2*b*a)/3,0, 10*exp(4*b*a)/9 - 2*exp(b*a)/9 + exp(-2*b*a)/9]])
    # 记录归一化因子
    max_val = np.max(np.abs(mat))
    # 仅当max_val不为0时归一化
    if max_val != 0:
        return mat / max_val, max_val
    return mat, 1

# 修改特征值乘积计算循环
for i, b in enumerate(b_s):
    for j, a in enumerate(a_s):
        mat, max_val = foo(a,b)
        w, v = np.linalg.eig(mat)
        # 还原行列式：归一化后的行列式 * (max_val^3)
        store_det_by_eigen[i,j] = (w[0]*w[1]*w[2]) * (max_val ** 3)

3. 使用更高精度的浮点数类型

如果平台支持，使用np.float128提升计算精度，减少溢出和舍入误差：

def foo(a,b):
    # 指定float128类型
    return np.array([[exp(-2*b*a), 0,-exp(b*a)/3 + exp(-2*b*a)/3],
                    [0, exp(-2*b*a), 0],
                    [-exp(b*a)/3 + exp(-2*b*a)/3,0, 10*exp(4*b*a)/9 - 2*exp(b*a)/9 + exp(-2*b*a)/9]], dtype=np.float128)

注意：float128在Windows等部分平台可能不支持，且计算速度会有所下降。

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