厄米协方差矩阵部分块奇异无法求逆的排查与解决

问题背景 #

我正在复现一篇天文学论文，需要将厄米协方差矩阵求逆得到Fisher矩阵。构建的协方差矩阵存储在形状为(n_z, n_k, n_mu, 3, 3)的Gamma[i,m,n]数组中，其中n_z、n_k、n_mu分别对应红移z、波数k、方向余弦μ维度的数据点数量，每个内部3×3块都是厄米矩阵且理论上可逆。

已验证所有矩阵块的厄米性（对角元实数、非对角元互为共轭），但对每个块求逆时，900个块中有23个抛出LinAlgError: Singular matrix错误，这些块的行列式在浮点精度下为0。非奇异块的特征值存在微小虚部和负实部（数值噪声），但奇异块的行列式严格为0。矩阵构建基于CAMB库，核心代码如下：

import camb
from camb import model
import numpy as np
from scipy.integrate import quad

# --- Constants and conversions ---
c_light = 2.998e5  # km/s

h0 = 0.6774
Om = 0.31
Ob = 0.05

H0 = h0 * 100
ns = 0.967
As = 2.142e-9

# --- Cosmology functions (H, chi, Vsur, Nk, Pkm) ---
# (definitions omitted here for brevity; full code used in my script)

# --- CAMB setup ---
pars = camb.CAMBparams()
pars.set_cosmology(H0=H0, ombh2=Ob*h0**2, omch2=(Om-Ob)*h0**2, omk=0, mnu=0)
pars.InitPower.set_params(ns=ns, r=0, As=As)
pars.set_matter_power(redshifts=[0.133, 0.3, 0.467], kmax=0.2)
pars.NonLinear = model.NonLinear_none
results = camb.get_results(pars)

# --- Arrays ---
S_area = 10000
omega = S_area*(np.pi/180)**2
z = np.array([0.133, 0.3, 0.467])
Dz = 0.111

deltak = [kmin(zi, Dz, omega) for zi in z]
k = [np.logspace(np.log10(dk), np.log10(0.2), num=30) for dk in deltak]
k = np.array(k)                  # -> shape (n_z, n_kpoints)
mu = np.array([np.linspace(-1, 1, num=10) for _ in z])

Deltamu = 2
n_z = 3
n_k = 30
n_mu = 10

pk = np.array([Pkm(ki, zi) for zi in z for ki in k])  # ensure pk[i,m] is scalar in my loop
# ... compute f, h, biases, alphas, n_g arrays ...

Gamma = np.zeros((n_z, n_k, n_mu, 3, 3), dtype=complex)

def P_auto_tilde(mui, hi, ki, alpha, b, fi, ng, pki):
    return ((b + fi*mui**2)**2 + (hi/ki)**2 * alpha**2 * fi**2 * mui**2) * pki + ng

def Pxy(mui, hi, ki, a1, a2, b1, b2, fi):
    return (
        (b1 + fi*mui**2)*(b2 + fi*mui**2)
        + (hi/ki)**2 * a1*a2 * fi**2 * mui**2
        - 1j * (fi*hi*mui*(a1*(b2+fi*mui**2) - a2*(b1+fi*mui**2)) / ki)
    )

for i in range(n_z):
    for m in range(n_k):
        for n in range(n_mu):
            mu_val = mu[i, n]
            h_val  = h[i]
            k_val  = k[i, m]
            f_val  = f[i]
            pk_val = pk[i, m]

            Pxx = P_auto_tilde(mu_val, h_val, k_val, a1_val, b1_val, f_val, ng1_val, pk_val)
            Pyy = P_auto_tilde(mu_val, h_val, k_val, a2_val, b2_val, f_val, ng2_val, pk_val)
            Pxy_val = Pxy(mu_val, h_val, k_val, a1_val, a2_val, b1_val, b2_val, f_val) * pk_val
            Pyx_val = Pxy_val.conj()

            pref = 2.0 / Nk(z[i], k_val, deltak[i], Deltamu, Dz, omega)

            M = np.zeros((3, 3), dtype=complex)
            M[0, 0] = Pxx**2
            M[0, 1] = Pxx * Pxy_val
            M[0, 2] = Pxy_val**2

            M[1, 0] = M[0, 1].conj()
            M[1, 1] = 0.5 * (Pxx*Pyy + Pxy_val*Pyx_val)
            M[1, 2] = Pxy_val * Pyy

            M[2, 0] = M[0, 2].conj()
            M[2, 1] = M[1, 2].conj()
            M[2, 2] = Pyy**2

            Gamma[i, m, n] = M * pref

奇异原因排查 #

1. 矩阵结构的代数相关性 #

观察3×3矩阵M的构造，其行/列存在潜在线性相关性：
当满足Pxx·Pyy = |Pxy|²时，矩阵的秩会降为1或2，直接导致奇异。可以对奇异块验证abs(Pxx*Pyy - Pxy_val * Pyx_val)是否在浮点精度内接近0，若成立则说明两个自功率谱与交叉功率谱完全相关，引发秩退化。

2. 数值计算精度问题 #

• 当Pxx或Pyy非常小时（比如小k值下物质功率谱P(k)趋近于0，且噪声项ng设置不合理），Pxx²或Pyy²会因浮点下溢变为0，导致矩阵对角元为0。
• 检查Nk函数的输出：若Nk数值极大，pref因子会让Gamma所有元素趋近于0，引发数值奇异。

3. 物理模型与采样的问题 #

• 当mu=0时，mui²=0，P_auto_tilde和Pxy的结构会简化，可能导致矩阵秩退化。检查奇异块是否集中在mu=0附近。
• 若参数a1=a2且b1=b2，矩阵会呈现高度对称结构，也可能引发秩退化。

解决方法（不使用pinv()的正则化方案） #

1. 物理意义正则化 #

• 添加自适应噪声项：在Gamma矩阵的对角元上添加极小的正实数ε，取值需远小于物理信号量级，比如ε = 1e-12 * np.mean(np.diag(Gamma[i,m,n]))，基于对角元均值缩放，确保矩阵正定且不影响物理结果。
• 修正功率谱计算：检查小k区域Pkm函数的实现，若P(k)过小，可使用CAMB的线性外推功能避免数值下溢。

2. 数值正则化技巧 #

• 特征值截断修正：利用厄米矩阵的特征值分解（np.linalg.eigh），将小于阈值的特征值替换为阈值，重构矩阵后求逆：

for i in range(n_z):
    for m in range(n_k):
        for n in range(n_mu):
            mat = Gamma[i,m,n]
            try:
                inv_mat = np.linalg.inv(mat)
            except np.linalg.LinAlgError:
                vals, vecs = np.linalg.eigh(mat)
                # 替换极小特征值，阈值可根据数据调整
                vals[vals < 1e-10] = 1e-10
                reg_mat = vecs @ np.diag(vals) @ vecs.conj().T
                inv_mat = np.linalg.inv(reg_mat)
            # 保存逆矩阵到结果数组

• 对角加载：在原矩阵上添加ε * I（I为3×3单位矩阵），简单直接且保证矩阵正定：

epsilon = 1e-12  # 根据实际数据量级调整
for i in range(n_z):
    for m in range(n_k):
        for n in range(n_mu):
            mat = Gamma[i,m,n]
            reg_mat = mat + epsilon * np.eye(3)
            inv_mat = np.linalg.inv(reg_mat)

3. 数据采样调整 #

检查奇异块对应的k、μ采样点：若集中在小k或μ=±1等物理退化区域，可调整采样密度，合并邻近bin或跳过这些区域，从源头上减少奇异情况。

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