Streaming Polyphase Factor Decomposition / SKAT

Download original file

This is the original source text. Reading it does not run or verify the program.

Read the typeset paper

1 import math
2 import unittest
3 import numpy as np
4 
5 from residual_trace_cap import (
6     residual_spectrum_certificate,
7     trace_cap_log_mgf_upper,
8     trace_frobenius_cap_log_mgf_upper,
9     trace_cap_tail_bracket,
10 )
11 
12 
13 class ResidualTraceCapTests(unittest.TestCase):
14     def test_exact_additive_identity(self):
15         rng = np.random.default_rng(1)
16         A = rng.normal(size=(18, 14))
17         U, _ = np.linalg.qr(rng.normal(size=(18, 7)), mode='reduced')
18         B = U.T @ A
19         F = A - U @ B
20         lhs = A.T @ A
21         rhs = B.T @ B + F.T @ F
22         rel = np.linalg.norm(lhs-rhs) / np.linalg.norm(lhs)
23         self.assertLess(rel, 1e-12)
24 
25     def test_trace_is_exact(self):
26         rng = np.random.default_rng(2)
27         A = rng.normal(size=(20, 15))
28         U = rng.normal(size=(20, 8))
29         cert = residual_spectrum_certificate(A, U)
30         true_trace = float(np.trace(A.T @ A))
31         self.assertAlmostEqual(cert.exact_total_trace, true_trace, places=10)
32 
33     def test_log_mgf_envelope_dominates_exact(self):
34         rng = np.random.default_rng(3)
35         A = rng.normal(size=(16, 12))
36         U = rng.normal(size=(16, 6))
37         cert = residual_spectrum_certificate(A, U)
38         H = A.T @ A
39         lam = np.linalg.eigvalsh(H)
40         theta = 0.1 / np.linalg.eigvalsh(H).max()
41         exact = float(-0.5*np.sum(np.log1p(-2*theta*lam)))
42         upper = trace_cap_log_mgf_upper(theta, cert)
43         self.assertLessEqual(exact, upper + 1e-10)
44 
45     def test_tail_bracket_contains_exact_kcos_tail(self):
46         from kcos_v2 import adaptive_kcos_logsf
47         rng = np.random.default_rng(4)
48         A = rng.normal(size=(18, 14))
49         U = rng.normal(size=(18, 9))
50         cert = residual_spectrum_certificate(A, U)
51         lam = np.linalg.eigvalsh(A.T @ A)[::-1]
52         mean = float(np.sum(lam))
53         q = mean + 5*math.sqrt(2*np.sum(lam*lam))
54         true = adaptive_kcos_logsf(q, lam, relative_tolerance=1e-9)
55         br = trace_cap_tail_bracket(q, cert, relative_tolerance=1e-8)
56         self.assertLessEqual(br.lower_log_p, true.log_p + 1e-9)
57         self.assertGreaterEqual(br.upper_log_p, true.log_p - 1e-9)
58 
59     def test_trace_frobenius_cap_envelope_dominates_exact(self):
60         rng = np.random.default_rng(5)
61         for _ in range(12):
62             A = rng.normal(size=(15, 11))
63             U = rng.normal(size=(15, 7))
64             cert = residual_spectrum_certificate(A, U)
65             H = A.T @ A
66             lam = np.linalg.eigvalsh(H)
67             theta = 0.08 / lam.max()
68             exact = float(-0.5*np.sum(np.log1p(-2*theta*lam)))
69             upper = trace_frobenius_cap_log_mgf_upper(theta, cert)
70             self.assertLessEqual(exact, upper + 1e-10)
71 
72 
73 if __name__ == '__main__':
74     unittest.main(verbosity=2)

Back to source controls

Back to Streaming Polyphase Factor Decomposition / SKAT