From 4e443f95b29233a03e9cbb496bd645e469973fed Mon Sep 17 00:00:00 2001 From: Nick Petalas Date: Tue, 3 Nov 2026 00:45:46 +0000 Subject: [PATCH] feat(research): prove the full Gaussian inverse rank-error law --- research/gaussian-inverse-rank-error/PROOF.md | 176 ++++++++++++ .../gaussian-inverse-rank-error/README.md | 140 ++++++++++ .../gaussian-inverse-rank-error/SOURCES.json | 24 ++ .../gaussian-inverse-rank-error/TASK.json | 93 +++++++ .../gaussian-inverse-rank-error/VERIFY.json | 36 +++ .../gaussian-inverse-rank-error/verify.py | 250 ++++++++++++++++++ 6 files changed, 719 insertions(+) create mode 100644 research/gaussian-inverse-rank-error/PROOF.md create mode 100644 research/gaussian-inverse-rank-error/README.md create mode 100644 research/gaussian-inverse-rank-error/SOURCES.json create mode 100644 research/gaussian-inverse-rank-error/TASK.json create mode 100644 research/gaussian-inverse-rank-error/VERIFY.json create mode 100644 research/gaussian-inverse-rank-error/verify.py diff --git a/research/gaussian-inverse-rank-error/PROOF.md b/research/gaussian-inverse-rank-error/PROOF.md new file mode 100644 index 0000000000..71b623a637 --- /dev/null +++ b/research/gaussian-inverse-rank-error/PROOF.md @@ -0,0 +1,176 @@ +# Exact fixed-rank compression error of the full Gaussian inverse + +This contribution was produced by Nick Petalas. Apache-2.0. + +## Claim and exact predecessor premise + +Fix alpha>0, put b=pi/alpha^2, and take 0 0) d_r(K_theta^(-1)) / ||K_theta^(-1)||_op + = 1 / sqrt(2 floor(r/2)+1). (1) + +In particular, a rank-one approximation has asymptotically the same +operator-norm error as the zero approximation. Rank two first lowers the +relative limiting error to 1/sqrt(3). This statement allows arbitrary +linear approximants, including those mixing the two ends of the cell. +It is an analytic optimum, not a method for preparing or applying modes. + +The cited premise is JoyboyBrian's [PR5183](https://github.com/CrocSwap/integer-mult-bounds/pull/5183), +retained at base 7c49fe32c55de02c2c873655c804b16e4c3fef0a in +[PROOF.md](../gaussian-half-line-operator-limit/PROOF.md). Its theorem gives +an isometric embedding U_theta into H=L2(0,infinity) direct-sum L2(0,infinity) +and the checkerboard sign matrix S=diag((-1)^i) such that + + X_theta = U_theta [sqrt(theta) S K_theta^(-1) S] U_theta* + ||X_theta - a (T_b direct-sum T_b)||_op -> 0, + a = alpha exp(b/4) / P(b)^3, + P(b) = product_(j>=1) (1-exp(-2bj)), + (T_b f)(x) = integral_(x to infinity) exp(-b(y^2-x^2)) f(y) dy. (2) + +We use (2) as a credited mathematical premise. We do not re-prove or +replay the source finite Gaussian matrices, source exponents, q-LDU sum, +embedding estimates or boundary-layer convergence. The new argument +calculates every singular value of T_b and transfers the resulting +approximation law through (2). No uniform growing-alpha result is used. + +## Complete half-line singular system + +Let H_n be the physicists' Hermite polynomial. Define the normalized +full-line Hermite functions, scaled to b, by + + phi_n(x) = (2b)^(1/4) H_n(sqrt(2b)x) exp(-bx^2) + / [pi^(1/4) sqrt(2^n n!)]. + +Their orthonormality, completeness and parity are classical; see +[NIST DLMF 1.18.42-43](https://dlmf.nist.gov/1.18#E42). +Let psi_n be sqrt(2) times the restriction of phi_n to x>0. +Each of the two sequences + + e_k = psi_(2k), o_k = psi_(2k+1), k>=0, + +is an orthonormal basis of the SAME half-line Hilbert space. To see +completeness, extend a half-line function evenly or oddly to the full +line, use the full-line basis, and restrict back. The two sequences are +not asserted to be mutually orthogonal to each other. + +The classical weighted derivative identity is + + d/dz [exp(-z^2) H_n(z)] = -exp(-z^2) H_(n+1)(z). + +See [NIST DLMF 18.9.26](https://dlmf.nist.gov/18.9#E26). +The polynomial times Gaussian vanishes at infinity. With z=sqrt(2b)y, +integrating this identity gives + + T_b [H_(2k+1)(sqrt(2b)y) exp(-by^2)](x) + = H_(2k)(sqrt(2b)x) exp(-bx^2) / sqrt(2b). + +Taking account of the different normalization constants yields + + T_b o_k = s_k e_k, + s_k = 1 / [2 sqrt(b(2k+1))]. (3) + +There is no extra boundary condition imposed on either input or output. +The integral itself determines (3). + +For completeness T_b is bounded independently of the spectral conclusion. +For 0<=x<=y one has y^2-x^2 >= (y-x)^2. Both integrals of its nonnegative +kernel, in either variable with the other fixed, are at most +integral_0^infinity exp(-bt^2)dt=sqrt(pi)/(2sqrt(b)). The Schur test +therefore bounds its operator norm. Equality (3) on the input basis +extends by continuity to every L2 input. Since the output sequence is +also a basis, this is a complete singular-value expansion. Its adjoint +maps e_k to s_k o_k. Truncating after m terms has norm error s_m -> 0, +so T_b is compact. Its singular values, indexed from zero, are exactly +the strictly decreasing numbers s_k in (3). + +The Hermite facts and diagonal-operator argument are classical. The new +source-specific consequence is the full spectrum and fixed-rank law for +(2), rather than another estimate of its top norm. + +## Optimal rank error and doubled modes + +Put B=a(T_b direct-sum T_b). Each s_k occurs twice in the singular +spectrum of B, multiplied by a. Thus + + d_r(B) = a / [2 sqrt(b(2 floor(r/2)+1))]. (4) + +Here is a direct minimax proof, including arbitrary approximants. The +upper bound retains the r largest singular modes. For the lower bound, +restrict any rank-at-most-r operator R to the span of the first r+1 +right singular vectors. Rank-nullity supplies a nonzero vector v in +that span with Rv=0. Orthogonality of their images implies +||Bv|| >= sigma_(r+1)(B)||v||, with singular values now indexed from one. +Consequently ||B-R|| >= sigma_(r+1)(B), which is (4). A rank mixing the +two edges cannot escape this argument. + +For any bounded X,Y on the same spaces, + + |d_r(X)-d_r(Y)| <= ||X-Y||. (5) + +Indeed, for every rank-at-most-r R, the triangle inequality gives +||X-R|| <= ||X-Y||+||Y-R||; take infima and exchange X,Y. +Applying (5) to (2) and (4) gives the absolute limit + + sqrt(theta) d_r(K_theta^(-1)) + -> a / [2 sqrt(b(2 floor(r/2)+1))]. (6) + +An isometric embedding followed by zero on its orthogonal complement +does not change nonzero singular values or rank approximation errors. +One may also prove this directly: embed a finite approximant for the +upper bound and compress any Hilbert-space approximant by U_theta* +and U_theta for the lower bound. Compression does not increase rank +or norm. Multiplication by S on either side preserves both. These facts +justify the left side of (6) without changing the permitted approximants. +For r=0, (6) has the positive limit a/(2sqrt(b)). Dividing proves (1). + +Equivalently, for every fixed j>=1, + + sqrt(theta) sigma_j(K_theta^(-1)) + -> a / [2 sqrt(b(2 floor((j-1)/2)+1))]. (7) + +The order of limits is fixed j or r first, then theta tending to zero. + +## Exact compression threshold in the limiting operator + +For 00 +and let theta tend to zero. For every fixed nonnegative integer r, + +```text +best rank-at-most-r approximation error / inverse operator norm + -> 1 / sqrt(2 floor(r/2)+1). +``` + +This is a scoped written mathematical proof with exact finite corroboration. +The new result determines the complete limiting singular spectrum and the +optimal fixed-rank error, beyond the predecessor's top-norm result. +The two boundary layers duplicate every singular value. Consequently rank +one cannot lower the limiting relative error, and the limiting operator +needs rank four to reach relative error at most one half. These are analytic +rank counts, not physical recursive bills. No new kappa is claimed. + +[PROOF.md](PROOF.md) defines the exact finite matrix and proves the result. +It constructs odd/even half-line Hermite bases and derives +`s_k(T_b)=1/(2 sqrt(b(2k+1)))`. Rank-nullity proves optimality even for +approximants mixing the two edges. The source operator-norm limit and the +Lipschitz property of rank approximation errors transfer the result to the +finite Gaussian inverse. The source theorem is a credited premise; no +predecessor scientific replay is included. + +## Reproduction + +From the repository root, with CPython 3.12 and its standard library: + +```sh +/Users/nick/git/integer-mult-bounds/.venv/bin/python3.12 -I -W error -B research/gaussian-inverse-rank-error/verify.py +/Users/nick/git/integer-mult-bounds/.venv/bin/python3.12 -I -W error -O -B research/gaussian-inverse-rank-error/verify.py +/Users/nick/git/integer-mult-bounds/.venv/bin/python3.12 -B scripts/research_package.py --package research/gaussian-inverse-rank-error +/Users/nick/git/integer-mult-bounds/.venv/bin/python3.12 -B scripts/check_research_commit.py --package research/gaussian-inverse-rank-error --check-source-history +``` + +On another host, substitute its Python 3.12 path. The final command requires +a commit containing this package and full source history. The direct commands +print `PASS_GAUSSIAN_INVERSE_RANK_ERROR` inside a generated JSON receipt bound +to the package and source-proof hashes. Optional `--receipt /absolute/path.json` +creates that receipt exclusively outside the checkout. No replay writes +artifacts into the checkout. + +The exact verifier compares explicit factorial and recurrence definitions of +Hermite polynomials, evaluates their Gaussian Gram integrals, and solves +weighted antiderivatives by descending coefficients. Coherent polynomial +mixtures test the integral map and normalization. Finite doubled diagonal +truncations check rank indexing, and rank-one matrices mixing the first two +edge modes have explicit kernel witnesses. Rational tolerance checks exercise +the least-rank threshold, including its equality cases. Invalid contracts and +one-edge, parity and dilation controls remain active under optimized Python. +These are finite corroboration, not a substitute for the written proof or +classical Hermite completeness. + +Local CPython 3.12.15 on macOS passed both registered commands. Individually +timed direct runs took 0.137 and 0.146 seconds, with identical JSON receipts. +They enumerate 34 polynomial definitions, 33 derivative identities, 512 +Gram entries, 16 antiderivatives and normalizations, 81 coherent mixtures, +1,275 coordinate truncations, 324 rational tolerances, 33 mixed-edge +rank-one matrices and 15 rejection/semantic controls. These are check +durations and finite counts, not a runtime speedup claim. + +## Scope and limits + +The order of limits is fixed r or fixed singular index first, then theta +shrinking at fixed alpha. The exact limiting-operator least rank is +`2 ceil((epsilon^(-2)-1)/2)` for `0 None: + if not condition: + raise ValueError(message) + + +def trim(p: tuple[F, ...]) -> tuple[F, ...]: + for j in range(len(p) - 1, -1, -1): + if p[j]: + return p[:j + 1] + return (F(0),) + + +def add(p: tuple[F, ...], q: tuple[F, ...]) -> tuple[F, ...]: + size = max(len(p), len(q)) + return trim(tuple((p[i] if i < len(p) else F(0)) + + (q[i] if i < len(q) else F(0)) for i in range(size))) + + +def scale(p: tuple[F, ...], c: F) -> tuple[F, ...]: + return trim(tuple(c * v for v in p)) + + +def derivative(p: tuple[F, ...]) -> tuple[F, ...]: + return trim(tuple(i * p[i] for i in range(1, len(p)))) + + +def hermite_explicit(n: int) -> tuple[F, ...]: + require(type(n) is int and 0 <= n <= MAX_DEGREE, 'invalid degree') + p = [F(0)] * (n + 1) + for j in range(n // 2 + 1): + p[n - 2 * j] = F((-1)**j * math.factorial(n) * 2**(n - 2 * j), + math.factorial(j) * math.factorial(n - 2 * j)) + return tuple(p) + + +def hermite_recurrence() -> tuple[tuple[F, ...], ...]: + terms = [(F(1),), (F(0), F(2))] + for n in range(1, MAX_DEGREE): + terms.append(add((F(0),) + scale(terms[-1], F(2)), + scale(terms[-2], F(-2 * n)))) + require(len(terms) == MAX_DEGREE + 1, 'recurrence capacity') + return tuple(terms) + + +def reverse_integral(p: tuple[F, ...]) -> tuple[F, ...]: + """Solve 2x Q-Q'=P by descending coefficients, not a Hermite recurrence.""" + require(1 < len(p) <= MAX_DEGREE + 1 and len(p) % 2 == 0, + 'odd degree required') + require(all(not p[i] for i in range(0, len(p), 2)), 'odd parity required') + residual = list(p) + q = [F(0)] * (len(p) - 1) + for j in range(len(p) - 2, -1, -2): + q[j] = residual[j + 1] / 2 + residual[j + 1] -= 2 * q[j] + if j: + residual[j - 1] += j * q[j] + require(not any(residual), 'antiderivative remainder') + return trim(tuple(q)) + + +def gram(p: tuple[F, ...], q: tuple[F, ...]) -> F: + """Integral_0^infinity P Q exp(-x^2) dx / sqrt(pi), same parity only.""" + result = F(0) + for i, a in enumerate(p): + for j, b in enumerate(q): + if a and b: + require((i + j) % 2 == 0, 'mixed parity integral') + m = (i + j) // 2 + moment = F(math.factorial(2 * m), 2 * 4**m * math.factorial(m)) + result += a * b * moment + return result + + +def spectrum_square(b: F, rank: int, copies: int = 2) -> F: + require(type(b) is F and 0 < b <= 16, 'invalid b') + require(type(rank) is int and 0 <= rank <= 4096, 'invalid rank') + require(type(copies) is int and 1 <= copies <= 2, 'invalid copies') + return 1 / (4 * b * (2 * (rank // copies) + 1)) + + +def threshold(epsilon: F) -> int: + require(type(epsilon) is F and F(1, 32) <= epsilon <= 1, 'invalid tolerance') + x = (1 / epsilon**2 - 1) / 2 + return 2 * ((x.numerator + x.denominator - 1) // x.denominator) + + +def polynomial_checks() -> dict[str, int]: + hs = hermite_recurrence() + for n, p in enumerate(hs): + require(p == hermite_explicit(n), 'independent Hermite formula') + for n in range(MAX_DEGREE): + weighted = add((F(0),) + scale(hs[n], F(2)), scale(derivative(hs[n]), F(-1))) + require(weighted == hs[n + 1], 'weighted derivative') + for k in range(MODES): + require(reverse_integral(hs[2*k + 1]) == hs[2*k], 'literal inverse integral') + for parity in range(2): + n = 2*k + parity + for j in range(MODES): + expected = F(2**n * math.factorial(n), 2) if j == k else F(0) + require(gram(hs[n], hs[2*j + parity]) == expected, 'Gaussian Gram') + ratio = gram(hs[2*k], hs[2*k]) / gram(hs[2*k+1], hs[2*k+1]) + require(ratio == spectrum_square(F(1, 2), k, 1), 'singular normalization') + return {'Hermite_polynomials': len(hs), 'weighted_derivatives': MAX_DEGREE, + 'literal_antiderivatives': MODES, 'Gaussian_Gram_entries': 2*MODES**2, + 'singular_normalizations': MODES} + + +def mixture_checks() -> int: + hs = tuple(hermite_explicit(n) for n in range(8)) + count = 0 + for cs in itertools.product((-1, 0, 1), repeat=4): + p, q = (F(0),), (F(0),) + for k, c in enumerate(cs): + p = add(p, scale(hs[2*k + 1], F(c))) + q = add(q, scale(hs[2*k], F(c))) + if any(cs): + require(reverse_integral(p) == q, 'coherent integral mixture') + predicted = sum(F(c*c * 2**(2*k) * math.factorial(2*k), 2) + for k, c in enumerate(cs)) + require(gram(q, q) == predicted, 'coherent output norm') + require(gram(q, q) <= gram(p, p) / 2, 'contraction on mixture') + else: + require(q == (F(0),), 'zero mixture') + count += 1 + return count + + +def rank_checks() -> dict[str, int]: + subset_count = 0 + for b in (F(1, 8), F(1, 4), F(1, 2), F(1), F(7, 3)): + weights = tuple(spectrum_square(b, r) for r in range(8)) + for r in range(8): + errors = [] + for retained in itertools.combinations(range(8), r): + errors.append(max(v for j, v in enumerate(weights) if j not in retained)) + subset_count += 1 + require(min(errors) == spectrum_square(b, r), 'best coordinate truncation') + tolerances = set(F(i, j) for j in range(1, 33) for i in range(1, j + 1)) + for e in sorted(tolerances): + rank = threshold(e) + require(spectrum_square(F(1, 4), rank) <= e**2, 'threshold sufficiency') + if rank: + require(spectrum_square(F(1, 4), rank - 1) > e**2, 'threshold necessity') + mixed = 0 + for a, b, c, d in itertools.product((-1, 0, 1), repeat=4): + if a*d == b*c: + if a or b: + v = (b, -a) + elif c or d: + v = (d, -c) + else: + v = (1, 0) + require(v != (0, 0) and a*v[0] + b*v[1] == 0 and + c*v[0] + d*v[1] == 0, 'rank-one mixed-edge kernel') + mixed += 1 + return {'coordinate_truncations': subset_count, 'exact_tolerances': len(tolerances), + 'mixed_edge_rank_one_matrices': mixed} + + +def controls() -> int: + bad = ((0, 0, 2), (F(0), 0, 2), (F(1), -1, 2), + (F(1), True, 2), (F(1), 4097, 2), (F(1), 1, 0)) + count = 0 + for args in bad: + try: + spectrum_square(*args) + except ValueError: + count += 1 + else: + raise ValueError('invalid spectrum accepted') + for e in (F(0), F(1, 33), F(2), 0.5): + try: + threshold(e) + except ValueError: + count += 1 + else: + raise ValueError('invalid tolerance accepted') + require(spectrum_square(F(1), 1) != spectrum_square(F(1), 1, 1), + 'one-edge mutant must differ') + require(spectrum_square(F(1, 2), 2, 1) != F(1, 2*(2+1)), 'parity mutant') + require(spectrum_square(F(1, 8), 2) != spectrum_square(F(1, 2), 2), + 'missing dilation mutant') + require(threshold(F(1, 2)) == 4 and threshold(F(1, 10)) == 100, + 'threshold fixtures') + require(spectrum_square(F(1, 4), 1) == 1, 'rank-one duplicated top') + return count + 5 + + +def file_hashes() -> dict[str, str]: + paths = sorted(PACKAGE.iterdir()) + require(1 <= len(paths) <= 16, 'package file count') + result = {} + for path in paths: + require(path.is_file() and not path.is_symlink(), 'unexpected package entry') + require(path.stat().st_size <= 1024*1024, 'package file capacity') + result[path.name] = hashlib.sha256(path.read_bytes()).hexdigest() + manifest = json.loads((PACKAGE / 'SOURCES.json').read_text()) + rows = [row for row in manifest['files'] if row['path'] == SOURCE] + require(len(rows) == 1, 'missing source theorem pin') + source = ROOT / SOURCE + require(source.stat().st_size <= 1024*1024, 'source file capacity') + digest = hashlib.sha256(source.read_bytes()).hexdigest() + require(digest == rows[0]['sha256'], 'source theorem bytes changed') + result[SOURCE] = digest + return result + + +def main() -> None: + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument('--receipt', type=Path) + args = parser.parse_args() + before = file_hashes() + counts = polynomial_checks() + counts['coherent_mixtures'] = mixture_checks() + counts.update(rank_checks()) + counts['rejection_and_semantic_controls'] = controls() + require(file_hashes() == before, 'proof inputs changed during verification') + receipt = {'status': 'PASS_GAUSSIAN_INVERSE_RANK_ERROR', 'counts': counts, + 'claim_type': 'scoped written proof with finite exact corroboration', + 'sha256': before} + encoded = json.dumps(receipt, sort_keys=True, indent=2) + '\n' + if args.receipt is not None: + path = args.receipt.resolve() + require(not path.is_relative_to(ROOT), 'receipt must be outside the checkout') + with path.open('x', encoding='utf-8') as stream: + require(stream.write(encoded) == len(encoded), 'incomplete receipt write') + print(encoded, end='') + + +if __name__ == '__main__': + main()