4ad3394Analytic entropy-rate prediction; parallelize pasts across processesJeremy Magland 1"""Analytic prediction of the entropy rate of z = round(h * x), x iid
2N(0, sigma^2), from the Fourier modes H(f) = sum_j h_j e^{-2 pi i f j}.
4Derivation, in three steps:
61. Szego / Kolmogorov. y = h * x is stationary Gaussian with power spectrum
7 S(f) = sigma^2 |H(f)|^2. Its one-step prediction error variance from the
8 exact past is the geometric mean of the spectrum,
10 sigma_inf^2 = exp( int_0^1 ln S(f) df ),
12 so in the high-resolution regime (sigma_inf >> 1 quantization step)
14 H_rate ~ 1/2 log2(2 pi e sigma_inf^2).
16 This fails in two ways tied to quantization: it ignores that only the
17 *quantized* past is observed, and it diverges to -inf wherever H(f) has
18 zeros or deep stopbands (the log integral blows up while the true entropy
19 stays >= 0).
212. Quantization noise floor. Model roundoff as additive dither: z_t ~ y_t +
22 u_t with u_t iid U(-1/2, 1/2), variance 1/12, independent of y. The best
23 linear one-step prediction of w = y + u has error variance
24 exp(int ln(S + 1/12)) (Szego again, on S_w = S + 1/12), and u_{t+1} being
25 unpredictable splits off exactly, leaving for y_{t+1} itself
27 s*^2 = exp( int_0^1 ln(S(f) + 1/12) df ) - 1/12.
29 The 1/12 inside the log is what keeps the integral finite at spectral
30 zeros; s*^2 >= sigma_inf^2 always (geometric means are superadditive),
31 and s*^2 -> 0 as sigma -> 0.
333. Coarse-quantization saturation. Given the quantized past, y_{t+1} ~
34 N(m, s*^2) with a conditional mean m that equidistributes modulo the
35 quantization grid (valid when the marginal spread of y is >> 1 bin).
36 The average discrete entropy of round(m + N(0, s^2)) over a uniform grid
37 offset equals *exactly* the differential entropy of the Gaussian-uniform
38 convolution (the standard dithered-quantization identity):
40 G(s) = h( N(0, s^2) + U(-1/2, 1/2) ) [bits],
42 which tends to 1/2 log2(2 pi e s^2) for s >> 1 and to 0 for s -> 0
43 instead of diverging negative.
45Final formula:
47 H_rate ~ G( sqrt( exp( int_0^1 ln(sigma^2 |H(f)|^2 + 1/12) df ) - 1/12 ) )
49Known approximations: the dither term is really deterministic roundoff, not
50independent noise; w is not Gaussian (linear prediction is not optimal); and
51for kernels whose conditional mean does not equidistribute mod 1 (e.g. the
52identity kernel, where it is constant) G slightly overestimates at small s.
53"""
55import math
57import numpy as np
58from scipy.special import ndtr
60TWO_PI_E = 2.0 * math.pi * math.e
63def spectrum(kernel, n=1 << 18):
64 """|H(f)|^2 on the rfft grid f = k/n, k = 0..n/2."""
65 return np.abs(np.fft.rfft(np.asarray(kernel, dtype=float), n)) ** 2
68def log_spectrum_mean(kernel, sigma, floor=1.0 / 12.0, n=1 << 18):
69 """Mean over f in [0, 1) of ln(sigma^2 |H(f)|^2 + floor).
71 Trapezoid on [0, 1/2] (the spectrum is symmetric). floor > 0 keeps the
72 integrand finite at zeros of H.
73 """
74 logS = np.log(sigma * sigma * spectrum(kernel, n) + floor)
75 return float((logS.sum() - 0.5 * (logS[0] + logS[-1])) / (logS.size - 1))
78def sigma_infinity(kernel, sigma):
79 """Szego one-step prediction error std of y = h * x from its exact past:
80 the root geometric mean of the spectrum, sigma * |c_lead| * prod max(1, |b_k|)
81 over the zeros b_k of the tap polynomial (exact, robust to zeros of H on
82 the unit circle, where the log integral is still convergent)."""
83 c = np.trim_zeros(np.asarray(kernel, dtype=float)[::-1], 'f')
84 b = np.roots(c) if c.size > 1 else np.array([])
85 return sigma * abs(c[0]) * float(np.prod(np.maximum(1.0, np.abs(b))))
88def gauss_uniform_entropy(s):
89 """G(s) = h(N(0, s^2) + U(-1/2, 1/2)) in bits — the exact average entropy
90 of round(c + N(0, s^2)) over a uniform grid offset c."""
91 s = float(s)
92 if s <= 0:
93 return 0.0
94 if s < 1e-3:
95 return _edge_constant() * s
96 dv = min(s / 8.0, 0.01)
97 v = np.arange(0.0, 0.5 + 8.0 * s + 1.0, dv)
98 g = ndtr((v + 0.5) / s) - ndtr((v - 0.5) / s)
99 term = np.where(g > 0, -g * np.log2(np.clip(g, 1e-300, None)), 0.0)
100 return float(2.0 * dv * (term.sum() - 0.5 * term[0]))
103_EDGE_C = None
106def _edge_constant():
107 """int_-inf^inf h2(Phi(t)) dt: small-s slope of G(s)."""
108 global _EDGE_C
109 if _EDGE_C is None:
110 t = np.linspace(-12.0, 12.0, 20001)
111 p = np.clip(ndtr(t), 1e-300, 1 - 1e-16)
112 h2 = -(p * np.log2(p) + (1 - p) * np.log2(1 - p))
113 _EDGE_C = float(np.trapezoid(h2, t))
114 return _EDGE_C
117def predict_entropy_rate(kernel, sigma):
118 """Both predictions of the entropy rate of z = round(h * x), in bits.
120 Returns a dict:
121 sigma_inf Szego prediction error std from the exact past
122 highres 1/2 log2(2 pi e sigma_inf^2) — valid when sigma_inf >~ 1
123 s_star prediction error std of y_next from the quantized past
124 corrected G(s_star) — stays valid at spectral zeros and coarse
125 quantization
126 """
127 s_inf = sigma_infinity(kernel, sigma)
128 highres = (0.5 * math.log2(TWO_PI_E * s_inf * s_inf)
129 if s_inf > 0 else -math.inf)
130 gm_w = math.exp(log_spectrum_mean(kernel, sigma))
131 s_star = math.sqrt(max(gm_w - 1.0 / 12.0, 0.0))
132 return {
133 'sigma_inf': s_inf,
134 'highres': highres,
135 's_star': s_star,
136 'corrected': gauss_uniform_entropy(s_star),
137 }