Почему моя эрмитова ковариационная матрица не обратима? ⇐ Python

[Anonymous]

Я следую статье, в которой используется эрмитова ковариационная матрица


и инвертирует ее для создания матрицы Фишера. Я строю свою ковариацию как Gamma[i,m,n], хранящуюся внутри большего массива формы (n_z, n_k, n_mu, 3, 3), где n_z, n_k и n_mu — это количество точек данных в z, k и mu. Каждый внутренний блок 3×3 должен быть эрмитовым и обратимым.
Я проверил эрмитичность с помощью:

Код: Выделить всё

for i in range(n_z):
for m in range(n_k):
for n in range(n_mu):
print(np.allclose(Gamma[i, m, n], Gamma[i, m, n].conj().T, atol=1e-10, rtol=0))

Это выведет True для всех индексов. Я также явно симметризировал:

Код: Выделить всё

A = np.zeros_like(Gamma)
for i in range(n_z):
for m in range(n_k):
for n in range(n_mu):
A[i, m, n] = 0.5*(Gamma[i, m, n] + Gamma[i, m, n].T.conj())
print(np.allclose(A[i, m, n], Gamma[i, m, n], atol=1e-10))

Опять правда для всех. Диагональные записи действительны (

Код: Выделить всё

+0j

) и недиагонали являются взаимно сопряженными, как и ожидалось.
Однако, когда я пытаюсь инвертировать каждый блок 3×3, я иногда получаю сообщение LinAlgError: Singular матрица. Пример вывода:

Код: Выделить всё

for l in range(n_z):
for m in range(n_k):
for n in range(n_mu):
g = Gamma[l, m, n]
try:
invg = np.linalg.inv(g)
except np.linalg.LinAlgError:
print(f"Matrix not invertible at (l={l}, m={m}, n={n})")
print("det:", np.linalg.det(g))

например, печатает:

Код: Выделить всё

Matrix not invertible at (l=1, m=3, n=7)
0j
Matrix not invertible at (l=1, m=5, n=3)
0j

Итак, для этих индексов det(g) == 0 (в пределах плавающей точки), следовательно, в единственном числе. Это происходит для 23 из 900 блоков (не для всех), что является неожиданным, поскольку каждый блок должен представлять собой действительную ковариационную (полу)положительно определенную матрицу.
Я проверял собственные значения с помощью:

Код: Выделить всё

for i in range(n_z):
for m in range(n_k):
for n in range(n_mu):
print(np.linalg.eigvals(Gamma[i, m, n]))

Например, в (0,0,0) я получаю:

Код: Выделить всё

[ 1.70674158e+09-3.5e-09j,  2.29744006e+01-8.1e-09j,
-1.01287448e-08+1.2e-08j ]

Маленькие мнимые части выглядят как числовой шум, а небольшое отрицательное действительное значение третьего собственного значения также соответствует округлению до нуля. Но определитель для (0,0,0):

Код: Выделить всё

(4783.194482059961+83.43640215104267j)

который имеет заметную мнимую часть. Однако такие определители по-прежнему позволяют мне инвертировать его, только когда я получаю определитель, равный 0, он становится необратимым.
Код, который строит Gamma (большинство констант исправлено; 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

Почему некоторые блоки 3×3 являются единственными (детерминант нулевой), хотя они выглядят эрмитовыми, и как я могу исправить или диагностировать эту проблему? Каковы надежные способы обеспечить регуляризацию обратимости?
Изменить: благодаря подсказке @NickODell кажется вероятным, что некоторые строки/столбцы численно линейно зависят от проблемных точек данных, поскольку они представляют собой «приблизительно линейное масштабирование». Как бы вы порекомендовали решить эту проблему без использования pinv()?

Подробнее здесь: https://stackoverflow.com/questions/797 ... invertible