diff --git a/tests/test_lanczos.py b/tests/test_lanczos.py index 952eaba..98e0fcd 100644 --- a/tests/test_lanczos.py +++ b/tests/test_lanczos.py @@ -62,6 +62,48 @@ def test_high_degree(): alpha, beta = lanczos(A, v0=v0, rtol=1e-7, deg=k, orth=0) assert np.all(~np.isnan(alpha)) and np.all(~np.isnan(beta)) +def test_degree_unbiased(): + from primate.diagonalize import lanczos + from scipy.linalg import eigvalsh_tridiagonal + np.random.seed(1234) + n = 150 + A = symmetric(n) + alpha, beta = np.zeros(n, dtype=np.float32), np.zeros(n, dtype=np.float32) + + ## In general not guaranteed, but with full re-orthogonalization it's likely (and deterministic w/ fixed seed) + ew = np.linalg.eigvalsh(A) + for k in range(5, 150, 5): + v0 = np.random.choice([-1.0, +1.0], size=A.shape[1]) + alpha, beta = lanczos(A, v0=v0, rtol=1e-7, deg=k, orth=15) + assert np.all(~np.isnan(alpha)) and np.all(~np.isnan(beta)) + rr = eigvalsh_tridiagonal(alpha, beta) + assert np.allclose(rr.mean(), ew.mean(), atol=0.30) + +def test_toeplitz(): + from scipy.linalg import toeplitz + from primate.diagonalize import lanczos + from scipy.linalg import eigvalsh_tridiagonal + np.random.seed(1234) + n = 150 + c = np.random.uniform(size=n, low=0, high=1.0) + A = toeplitz(c) + alpha, beta = np.zeros(n, dtype=np.float32), np.zeros(n, dtype=np.float32) + + ## In general not guaranteed, but with full re-orthogonalization it's likely (and deterministic w/ fixed seed) + ew = np.linalg.eigvalsh(A) + for k in range(25): + v0 = np.random.choice([-1.0, +1.0], size=A.shape[1]) + alpha, beta = lanczos(A, v0=v0, rtol=1e-7, deg=150, orth=15) + assert np.all(~np.isnan(alpha)) and np.all(~np.isnan(beta)) + rr = eigvalsh_tridiagonal(alpha, beta) + assert np.allclose(rr.mean(), ew.mean(), atol=np.ptp(ew)*0.05) + +## TODO: go back and re-test fttr on toeplitz +# def test_symmetric(): +# assert np.allclose(ew_sym, np.flip(-ew_sym)) + + + def test_quadrature(): from primate.diagonalize import lanczos, _lanczos from sanity import lanczos_quadrature @@ -86,11 +128,51 @@ def test_quadrature_methods(): assert np.allclose(nw_test1, nw_test2) from primate.diagonalize import lanczos, _lanczos + np.random.seed(1234) n = 50 A = symmetric(n) + v0 = np.random.uniform(size=A.shape[1]) + a, b = lanczos(A, v0=v0, deg=n) + a, b = a, np.append([0], b) + quad1 = np.sum(_lanczos.quadrature(a, b, n, 0).prod(axis=1)) + quad2 = np.sum(_lanczos.quadrature(a, b, n, 1).prod(axis=1)) + assert np.isclose(quad1, quad2, atol=quad1*0.05) + +def test_quadrature_toeplitz(): + from primate.diagonalize import lanczos, _lanczos + from scipy.linalg import toeplitz + np.random.seed(1234) + + A = toeplitz(np.arange(1,8)).astype(np.float32) + v0 = np.random.uniform(size=A.shape[1]) + n = A.shape[1] a, b = lanczos(A, deg=n) a, b = a, np.append([0], b) quad1 = np.sum(_lanczos.quadrature(a, b, n, 0).prod(axis=1)) quad2 = np.sum(_lanczos.quadrature(a, b, n, 1).prod(axis=1)) - assert np.isclose(quad1, quad2, atol=0.5) + assert np.isclose(quad1, quad2, atol=np.abs(quad1)*0.50) + ## ground truth + ew, ev = np.linalg.eigh(A) + mu_0 = np.sum(np.abs(ew)) + mu_0 * np.ravel(ev[0,:])**2 + + np.array([orth_poly(ew[0], i, mu_0, a, b)**2 for i in range(n)]) + _lanczos.quadrature(a, b, n, 0)[:,1] + + # np.linalg.norm([orth_poly(ew[0], i, mu_0, a, b) for i in range(n)]) + ## + # + +# ## to mimick +# arr = np.array([orth_poly(ew[0], i, mu_0, a, b)**2 for i in range(n)]) + + n = 50 + c = np.random.uniform(size=n, low=0, high=1) + A = toeplitz(c) + v0 = np.random.uniform(size=A.shape[1]) + a, b = lanczos(A, deg=n) + a, b = a, np.append([0], b) + quad1 = np.sum(_lanczos.quadrature(a, b, n, 0).prod(axis=1)) + quad2 = np.sum(_lanczos.quadrature(a, b, n, 1).prod(axis=1)) + assert np.isclose(quad1, quad2, atol=quad1*0.50)