如何提升QuTiP函数数值精度？解决相对熵计算发散问题

核心思路 #

问题根源是极小特征值被误判为0导致对数运算发散，而相对熵的差值有限说明这些极小值的贡献在差值中抵消了，可从特征值处理、自定义计算逻辑、高精度浮点运算三个方向入手解决。

具体解决方案 #

1. 自定义相对熵计算（手动过滤极小特征值） #

QuTiP的entropy_relative底层依赖特征值分解，你可以手动提取特征值，过滤掉接近0的极小值后再计算：

import qutip as qt
import numpy as np

def entropy_relative_high_prec(state1, state2, eps=1e-15):
    # 提取算子的实特征值
    ev1 = np.real(state1.eigenenergies())
    ev2 = np.real(state2.eigenenergies())
    
    # 过滤小于阈值的特征值（视为数值意义上的0）
    valid_mask = (ev1 > eps) & (ev2 > eps)
    ev1_filtered = ev1[valid_mask]
    ev2_filtered = ev2[valid_mask]
    
    # 按相对熵定义计算：sum(p_i * (ln p_i - ln q_i))
    return np.sum(ev1_filtered * (np.log(ev1_filtered) - np.log(ev2_filtered)))

eps可根据精度需求调整，比如用1e-16适配双精度浮点的极限精度。

2. 切换到高精度浮点运算 #

直接用NumPy的四精度浮点类型重新计算特征值，提升数值区分度：

import qutip as qt
import numpy as np
from scipy.linalg import eigvalsh

def entropy_relative_high_prec(state1, state2):
    # 转换算子矩阵为四精度浮点
    rho1 = state1.full().astype(np.float128)
    rho2 = state2.full().astype(np.float128)
    
    # 高精度特征值分解
    ev1 = eigvalsh(rho1)
    ev2 = eigvalsh(rho2)
    
    # 仅保留大于四精度eps的特征值
    eps = np.finfo(np.float128).eps
    ev1 = ev1[ev1 > eps]
    ev2 = ev2[ev2 > eps]
    
    # 计算相对熵
    return np.sum(ev1 * (np.log(ev1) - np.log(ev2)))

注意：四精度运算会增加计算耗时，适合小规模量子系统。

3. 正则化算子规避零特征值 #

给算子加入极小的单位矩阵正则项，避免出现严格零特征值，再用QuTiP的算子对数函数计算：

import qutip as qt

def entropy_relative_regularized(state1, state2, eps=1e-15):
    # 添加正则项并归一化
    rho1_reg = state1 + eps * qt.qeye(state1.dims[0])
    rho2_reg = state2 + eps * qt.qeye(state2.dims[0])
    rho1_reg = rho1_reg / rho1_reg.tr()
    rho2_reg = rho2_reg / rho2_reg.tr()
    
    # 按相对熵的迹定义计算
    return qt.tr(rho1_reg * (rho1_reg.logm() - rho2_reg.logm()))

正则化的eps要足够小，避免对原始结果产生可察觉的偏差。

验证方式 #

用自定义函数分别计算两个相对熵后取差值，对比原方法的结果，确认差值的一致性和数值稳定性。

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