3. 冯·诺依曼熵的计算与数值稳定性

冯·诺依曼熵是量子信息论中最基础的量之一，它把经典香农熵推广到量子态上。对一个密度矩阵 ρ，其定义为 S(ρ) = −Tr(ρ log₂ ρ)。这个式子看起来简洁，但真正落到工程实现时，会立刻遇到两个问题：一是 log 作用于矩阵而不是标量，二是当 ρ 存在接近零或严格为零的本征值时，log 会发散。本节围绕定义、性质、本征值分解、对数处理、小本征值截断与数值误差控制展开，目标是让读者既能写出正确的计算代码，也能判断结果是否可信。

定义与基本性质

密度矩阵 ρ 是半正定、迹为 1 的厄米矩阵。由于 ρ 可对角化，设其本征值为 λᵢ，则 S(ρ) = −Σᵢ λᵢ log₂ λᵢ。这里约定 0 log 0 = 0，这是熵定义中的标准极限约定，不是近似。冯·诺依曼熵具有几个关键性质：非负性 S(ρ) ≥ 0；纯态熵为 0；在 d 维系统中最大值为 log₂ d，当且仅当 ρ = I/d 时取到；对任意酉变换 U 有 S(UρU†) = S(ρ)；对复合系统的约化态，熵可以大于任一子系统的最大熵，这正是纠缠的体现。

这些性质在数值实现中非常有用。例如，非负性可以作为结果合理性的第一道检查；最大值 log₂ d 可以作为上界，一旦计算结果超过它，说明代码或输入一定有问题。纯态熵为 0 则提醒我们，当所有本征值中只有一个为 1、其余为 0 时，截断策略必须保证结果严格为 0，而不是因为浮点误差得到 −0 或极小负数。

本征值分解：计算的核心路径

直接对矩阵做 log 运算在数值上既慢又不稳定，工程上几乎总是走本征值分解这条路。对厄米矩阵 ρ，使用 eigh 而非通用 eig，因为 eigh 利用对称性，返回实本征值并保证正交归一的本征向量，数值精度更高。计算流程可以概括为：输入 ρ → 检查厄米性与迹 → 本征分解 → 得到 λᵢ → 对每个 λᵢ 计算 −λᵢ log₂ λᵢ → 求和。

下面是一段可直接运行的 Python 示例，使用 NumPy 实现，并加入了基本的数值保护：

import numpy as np

def von_neumann_entropy(rho, tol=1e-12):
    rho = np.asarray(rho, dtype=complex)
    # 1. 厄米性检查
    if not np.allclose(rho, rho.conj().T, atol=1e-10):
        raise ValueError("密度矩阵必须厄米")
    # 2. 迹检查
    tr = np.trace(rho).real
    if not np.isclose(tr, 1.0, atol=1e-10):
        raise ValueError(f"迹应为1，实际为{tr}")
    # 3. 本征分解（只取实部，eigh保证实本征值）
    eigvals = np.linalg.eigvalsh(rho)
    # 4. 小本征值截断：小于tol的视为0
    eigvals = eigvals[eigvals > tol]
    # 5. 归一化修正（截断后重新归一，控制误差）
    eigvals = eigvals / eigvals.sum()
    # 6. 熵计算
    return -np.sum(eigvals * np.log2(eigvals))

这段代码里，第 4 步和第 5 步是数值稳定性的关键。eigh 返回的本征值可能包含 −1e−16 这样的微小负数，直接取 log 会得到 NaN。截断阈值 tol 的选择需要结合具体问题的精度要求，常见取值在 1e−12 到 1e−15 之间。截断后重新归一化，是为了补偿被丢弃的微小概率，使剩余本征值之和仍为 1，避免熵被系统性低估。

对数处理与截断策略

对数处理的核心原则是：只对严格大于零的本征值取对数。在浮点运算中，本征值可能因为舍入误差变成极小的负数，例如 −2.3e−16。如果直接计算 log₂(λ)，会得到 NaN 或复数，整个熵值失效。因此，在取对数之前，必须先把所有小于阈值 tol 的本征值置零或丢弃。这个操作称为小本征值截断。

截断阈值 tol 的选取需要权衡。tol 太大，会丢弃真实存在的小概率分量，导致熵被低估；tol 太小，则无法滤除浮点噪声，可能引入 NaN。一个实用的经验是：tol 取机器精度 eps 的若干倍，例如 1e−14 到 1e−12。对于维度 d 较大的系统，本征值分布更密集，tol 可以适当放宽到 1e−10，但需要同时检查截断前后熵的变化量是否在可接受范围内。

截断之后，剩余本征值之和会略小于 1。此时应做一次重新归一化：λᵢ ← λᵢ / Σλⱼ。这一步的代价很小，但能显著改善熵的精度。如果不做归一化，当被截断的总概率达到 1e−6 时，熵的偏差可能达到 1e−5 量级，对于高精度要求的场景不可忽略。

数值误差控制与结果判读

除了截断和归一化，还有几个值得养成的工程习惯。第一，始终使用 eigh 而不是 eig，前者对本征值排序并保证实数输出。第二，在计算前检查 ρ 的厄米性和迹，很多错误其实来自输入矩阵本身不合法。第三，计算完成后做上界检查：如果 S 大于 log₂ d 加上一个小的容差，说明计算过程有问题；如果 S 是负数且绝对值超过容差，同样说明截断或归一化有误。第四，对于纯态，S 应当严格为 0，可以用这个特例做单元测试。

下面给出一个简单的验证示例，构造一个 2 维最大混合态和一个纯态，观察输出是否符合预期：

# 最大混合态：S = log2(2) = 1
rho_mixed = np.eye(2) / 2
print(von_neumann_entropy(rho_mixed))  # 应输出约 1.0

# 纯态：S = 0
psi = np.array([1, 0], dtype=complex)
rho_pure = np.outer(psi, psi.conj())
print(von_neumann_entropy(rho_pure))   # 应输出约 0.0

在实际工程中，还经常遇到本征值谱跨越多个数量级的情况，例如纠缠谱中前几个本征值在 1e−1 量级，尾部延伸到 1e−15。此时固定阈值可能不够灵活，可以考虑自适应阈值：先按降序排列本征值，找到第一个小于 eps 的位置，只保留其之前的分量。另一种做法是使用对数数值稳定的求和方式，例如先计算每个 −λᵢ log₂ λᵢ 项，再用 math.fsum 或 Kahan 求和减少累加误差。对于维度达到数百以上的系统，这些细节会直接影响熵的第五位有效数字。

总结来说，冯·诺依曼熵的计算本身并不复杂，难点在于数值稳定性。本征值分解提供了正确的数学路径，小本征值截断和对数保护解决了 log 发散问题，重新归一化补偿了截断带来的概率损失，而上界检查和特例测试则构成了结果可信度的最后一道防线。把这几个环节都处理好，就能在工程实现中稳定地得到可信的熵值。