3. 冯·诺依曼熵的计算与数值实现

冯·诺依曼熵是量子信息中最核心的度量之一,它把经典香农熵推广到量子态上,用来刻画一个量子态所携带的“不确定程度”或“混合程度”。本节从定义出发,推导它的基本性质,然后落到可执行的数值实现上,重点处理零本征值与数值精度这两个在实际计算中几乎必然遇到的问题。

一、定义与推导

设量子系统处于密度矩阵 ρ 描述的状态。ρ 是半正定、迹为 1 的厄米矩阵。冯·诺依曼熵定义为

S(ρ) = −Tr(ρ ln ρ)

由于 ρ 厄米,可以对角化为 ρ = Σᵢ λᵢ |i⟩⟨i|,其中 λᵢ 是本征值,满足 λᵢ ≥ 0 且 Σᵢ λᵢ = 1。代入定义,利用迹在正交基下的展开,得到

S(ρ) = −Σᵢ λᵢ ln λᵢ

这正是香农熵的形式,只不过“概率分布”换成了密度矩阵的本征值谱。这个等价关系是数值实现的出发点:计算熵不需要真的做矩阵对数,只需要求本征值再套用标量公式。对纯态,谱为 (1, 0, 0, …),熵为 0;对 d 维最大混合态 ρ = I/d,谱为 (1/d, …, 1/d),熵为 ln d。这两个端点给出了熵的取值范围 0 ≤ S ≤ ln d(以自然对数为单位,若用比特则除以 ln 2)。

几个常用性质值得记住:熵在酉演化下不变,因为酉变换不改变本征值谱;熵是凹函数,即 S(Σ pᵢρᵢ) ≥ Σ pᵢS(ρᵢ),混合会增大不确定性;对复合系统,S(ρ_AB) ≤ S(ρ_A) + S(ρ_B),即次可加性。这些性质在后续讨论纠缠熵时会反复用到,尤其是纯态二分系统的纠缠熵等于任一约化态的冯·诺依曼熵这一结论。

二、数值实现的关键细节

直接按公式写代码会遇到两个陷阱。第一,当某个 λᵢ = 0 时,λᵢ ln λᵢ 在数学上取极限为 0,但数值上 ln 0 = −∞,0 × (−∞) 会得到 NaN。第二,密度矩阵经过数值运算(如部分迹、矩阵乘法、浮点累加)后,本征值可能出现极小的负值,例如 −1e−16,这是舍入误差而非物理真实。若不加处理,ln 负数同样产生 NaN。

处理策略是明确的:对 λᵢ 小于某个阈值(如 1e−12)的项直接跳过,视为 0 贡献;对略大于阈值但极小的本征值正常计算,因为它们的贡献虽小却真实。同时,在求本征值前可对 ρ 做一次厄米化 (ρ + ρ†)/2,抑制非厄米误差;求完后把本征值裁剪到 [0, 1] 区间并重新归一化,保证谱的物理合法性。下面给出一个可直接使用的 Python 实现:

import numpy as np

def von_neumann_entropy(rho, tol=1e-12, base=np.e):
    rho = np.asarray(rho, dtype=complex)
    rho = (rho + rho.conj().T) / 2.0        # 厄米化
    eigvals = np.linalg.eigvalsh(rho)        # 厄米矩阵用 eigvalsh,返回实本征值
    eigvals = np.clip(eigvals, 0.0, None)    # 裁掉微小负值
    s = eigvals.sum()
    if s > 0:
        eigvals = eigvals / s                # 重新归一化
    mask = eigvals > tol
    lam = eigvals[mask]
    S = -np.sum(lam * np.log(lam))
    if base != np.e:
        S = S / np.log(base)
    return float(S)

这里用 eigvalsh 而不是 eigvals,是因为密度矩阵厄米,前者更快且保证返回实数本征值,避免复数比较带来的麻烦。阈值 tol 的选取需要权衡:太小则微小负值漏网,太大则把真实的小本征值误判为 0。对双精度、维度在几十以内的典型计算,1e−12 是稳妥的默认值;若矩阵维度很大或经过多次迭代,可适当放宽到 1e−10。

三、用示例验证实现

用几个解析结果已知的态做校验,是确认实现正确的最快方式。考虑一个两能级态 ρ = [[0.7, 0.2], [0.2, 0.3]],其本征值为 0.5 ± √(0.04 + 0.04) = 0.5 ± 0.2828,即约 0.7828 与 0.2172。代入公式可得 S ≈ 0.5327(自然对数单位)。最大混合态 I/2 应给出 ln 2 ≈ 0.6931,纯态 |0⟩⟨0| 应给出 0。若实现返回这些值,说明本征值求解与零值处理都正确。

更贴近课程主题的校验是纠缠熵。取 Bell 态 |Φ⁺⟩ = (|00⟩ + |11⟩)/√2,对第二个量子比特求部分迹得到约化密度矩阵 I/2,其冯·诺依曼熵应为 ln 2。这个例子同时检验了部分迹、本征值求解和熵计算三个环节,是端到端验证的好用例。若把 Bell 态换成乘积态 |00⟩,约化态为纯态,熵应为 0,两者对比能直观看出纠缠与熵的对应关系。

冯·诺依曼熵数值计算流程与零本征值处理 输入密度矩阵 ρ 厄米 · 半正定 · Tr=1 厄米化 (ρ + ρ†)/2 求本征值谱 eigvalsh(ρ) → λᵢ 裁剪与归一 λᵢ ← max(λᵢ,0),再归一 λᵢ > tol ? tol = 1e−12 计入 −λᵢ ln λᵢ 真实小本征值保留贡献 是 跳过该项 视为 0 贡献,避免 NaN 否 求和得 S(ρ) S = −Σ λᵢ ln λᵢ 解析校验点 纯态→0 · I/2→ln2 · Bell→ln2 关键:本征值谱决定熵;零值与微小负值必须显式处理,否则出现 NaN

四、精度与稳定性实践

除了零本征值,还有几处容易出问题的地方。其一是本征值求解器对简并谱的敏感性:当多个本征值接近相等时,数值误差可能让它们略微偏离,但熵对谱的连续依赖保证了总误差仍然很小,一般无需特殊处理。其二是维度较大时的计算成本:eigvalsh 的复杂度约为 O(d³),对 d 达到几百的约化密度矩阵仍可接受,但若在参数扫描中反复调用,应尽量复用中间结果,避免重复对角化。

其三是单位约定。物理文献常用自然对数,信息论常用以 2 为底,代码中应显式指定 base 参数并在结果旁标注单位,否则跨模块比较时容易差一个 ln 2 的因子。其四是负熵的防护:理论上 S ≥ 0,但数值上若归一化不严,可能出现 −1e−16 量级的负值,输出前可做一次 max(S, 0) 的截断,并在调试阶段记录被截断的次数,作为数值健康度的指标。

建议在实现中保留一个“诊断模式”:返回熵的同时输出最小本征值、被跳过的本征值个数和归一化前后的迹。这样在后续计算纠缠熵、互信息等复合量时,一旦结果异常,可以快速定位是物理态本身有问题,还是数值处理环节出了偏差。

把上述实现封装好后,它就成为后续所有关联度量的基础构件:纠缠熵是约化态的冯·诺依曼熵,互信息是 S(ρ_A) + S(ρ_B) − S(ρ_AB),相对熵是 Tr(ρ ln ρ − ρ ln σ)。这些量都依赖同一个稳定的熵函数,因此在本节把零值处理、厄米化和归一化做扎实,后面章节的数值实验才能可靠复现。