Add test for covariance eigenvalue equivalence

Author: fonnesbeck

pymc/gp/cov.py

@@ -651,8 +651,12 @@ def power_spectral_density(self, omega: TensorLike) -> TensorVariable:
 
         z = pt.sqrt(2 * alpha) * pt.sqrt(pt.dot(pt.square(omega), pt.square(ls)))
         coeff = 2.0 * pt.power(2.0 * np.pi * alpha, D / 2.0) * pt.prod(ls) / pt.gamma(alpha)
-        term_z = pt.power(z / 2.0, nu) * pt.kv(nu, z)
-
+        
+        # Handle singularity at z=0
+        safe_z = pt.switch(pt.eq(z, 0), 1.0, z)
+        term_z = pt.power(safe_z / 2.0, nu) * pt.kv(nu, safe_z)
+        term_z = pt.switch(pt.eq(z, 0), pt.gamma(nu) / 2.0, term_z)
+        
         return coeff * term_z
 
 

tests/gp/test_cov.py

@@ -17,6 +17,7 @@
 import pytensor
 import pytensor.tensor as pt
 import pytest
+import scipy.linalg
 
 from scipy.special import gamma, iv, kv
 
@@ -553,6 +554,36 @@ def test_psd(self):
         )
         npt.assert_allclose(true_1d_psd, test_1d_psd, atol=1e-5)
 
+    def test_psd_eigenvalues(self):
+        # Test PSD implementation using Szegő’s Theorem
+        alpha = 1.5
+        ls = 0.5
+        N = 500
+        L = 50.0
+        dx = L / N
+        X = np.linspace(0, L, N)[:, None]
+
+        with pm.Model():
+            cov = pm.gp.cov.RatQuad(1, alpha=alpha, ls=ls)
+
+        K = cov(X).eval()
+
+        evals = scipy.linalg.eigvalsh(K)
+        evals = np.sort(evals)[::-1]
+
+        freqs = np.fft.fftfreq(N, d=dx)
+        omegas = 2 * np.pi * freqs
+
+        psd = cov.power_spectral_density(omegas[:, None]).eval()
+
+        psd_scaled = psd.flatten() / dx
+        psd_sorted = np.sort(psd_scaled)[::-1]
+
+        rel_err = np.abs((evals - psd_sorted) / psd_sorted)
+        med_rel_err = np.median(rel_err)
+
+        assert med_rel_err < 0.1
+
 
 class TestExponential:
     def test_1d(self):