From 5193f368a77a1533abcf8e8bafc97486f85fcb02 Mon Sep 17 00:00:00 2001 From: pinterbanget Date: Thu, 22 Oct 2026 13:23:39 +0700 Subject: [PATCH] =?UTF-8?q?Independently=20check=20#4554's=20native=20band?= =?UTF-8?q?ed=20inverse:=20exact=20elimination=20bounds,=20Gaussian=20tail?= =?UTF-8?q?,=20exact=20periodic=20Woodbury=20solves=20(no=20new=20=CE=BA)?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-Authored-By: Claude Opus 5.5 --- .../.gitattributes | 1 + .../PR4554-NATIVE-BANDED-INVERSE.md | 116 +++++++ .../pr4554-banded-inverse-check/README.md | 57 ++++ .../pr4554-banded-inverse-check/SOURCES.json | 18 + .../pr4554-banded-inverse-check/TASK.json | 100 ++++++ .../pr4554-banded-inverse-check/VERIFY.json | 18 + .../pr4554-banded-inverse-check/verify.py | 316 ++++++++++++++++++ 7 files changed, 626 insertions(+) create mode 100644 research/pr4554-banded-inverse-check/.gitattributes create mode 100644 research/pr4554-banded-inverse-check/PR4554-NATIVE-BANDED-INVERSE.md create mode 100644 research/pr4554-banded-inverse-check/README.md create mode 100644 research/pr4554-banded-inverse-check/SOURCES.json create mode 100644 research/pr4554-banded-inverse-check/TASK.json create mode 100644 research/pr4554-banded-inverse-check/VERIFY.json create mode 100644 research/pr4554-banded-inverse-check/verify.py diff --git a/research/pr4554-banded-inverse-check/.gitattributes b/research/pr4554-banded-inverse-check/.gitattributes new file mode 100644 index 0000000000..e80440330f --- /dev/null +++ b/research/pr4554-banded-inverse-check/.gitattributes @@ -0,0 +1 @@ +PR4554-NATIVE-BANDED-INVERSE.md -text diff --git a/research/pr4554-banded-inverse-check/PR4554-NATIVE-BANDED-INVERSE.md b/research/pr4554-banded-inverse-check/PR4554-NATIVE-BANDED-INVERSE.md new file mode 100644 index 0000000000..8be8f9922d --- /dev/null +++ b/research/pr4554-banded-inverse-check/PR4554-NATIVE-BANDED-INVERSE.md @@ -0,0 +1,116 @@ +# A native bounded-word inverse for a near-identity band matrix + +Let A be an s by s complex matrix with bandwidth h: a_ij=0 for |i-j|>h. Suppose A=I+E and ||E||_infinity<=rho<=1/2. All norms below use complex absolute values. The input entries are explicitly stored bounded-word fixed-point numbers; this is not an assumed inverse oracle. Put H=h+1. Gaussian elimination without pivoting, followed by forward and backward substitution, computes an approximation to A^-1 b in O(s H^2) bounded-word operations, O(s H) space and an explicit absolute error at most + + 1000 H^4 u (||b||_infinity+1), + +provided local complex rounding errors are at most u<=1/(1000 H^2), and enough integral bits are reserved as below. Constants are deliberately loose. To obtain output error <=2^-P it suffices to choose + + u <= 2^-P / [1000 H^4 (||b||_infinity+1)]. + +Thus P+log2(||b||+1)+4 log2 H+O(1) fractional bits suffice. If P, log s, log(||b||+1)=O(log n) and h<=n^C for fixed C, all operations use a fixed number of O(log n)-bit words. + +## Exact elimination has bounded growth + +Initially every row has diagonal-dominance margin + + |a_ii|-sum_(j!=i)|a_ij| >= 1-rho >=1/2, + +and absolute row sum at most 1+rho<=3/2. Consider elimination of column k in a remaining row i. Write l=a_ik/a_kk. For the Schur row after deleting column k, the triangle inequality gives + + new margin >= old margin + |a_ik|[1-offsum(row k)/|a_kk|] + >= old margin. + +Similarly its absolute row sum is at most + + old row sum - |a_ik| + |l| offsum(row k) + <= old row sum. + +Consequently all exact pivots are at least 1/2, every surviving row sum is at most 3/2, and every multiplier has magnitude at most 3. These assertions hold for every leading principal subproblem too. The inverse infinity norm is at most 2: if |x_i|=||x||_infinity, then |(Ax)_i|>= (|a_ii|-sumoff)|x_i| >= ||x||/2. + +No pivoting is needed. The lower and upper bandwidths remain h: a row only interacts with at most h preceding pivots and at most h following columns. This proves the claimed operation and space counts directly. + +## Rounded elimination retains these bounds + +Round a multiplier within absolute error u, and each Schur entry update within absolute error u after its exact bounded-word product and subtraction. As long as remaining row sums are at most 2, replacing an exact multiplier by the rounded multiplier perturbs an updated entry by at most 2u, and its final rounding adds at most u. There are at most 2H active entries per updated row and at most H elimination steps affecting that row. The total adverse change in either its row sum or its diagonal-dominance margin is therefore at most 6 H^2 u. The assumption u<=1/(1000 H^2) makes this less than 1/100. By induction the row sums stay below 2, the margins stay above 1/4, pivots have magnitude at least 1/4 and multipliers are bounded by 8. This closes the estimates used in the induction. + +Denote the computed factors by Lhat,Uhat. The factor residual is small without exponentially multiplying rounding errors. For one fixed row i, expand its successive updates literally: + + current a_ij = original a_ij - sum_(k already eliminated) lhat_ik uhat_kj + + sum(local update roundings). + +Each entry has at most H local update roundings. At the removed pivot column, rounding its multiplier contributes at most |uhat_kk|u<=2u. Thus each entry of Lhat Uhat-A has magnitude at most (H+2)u, with at most 2H+1 nonzero entries per row. In particular + + ||Lhat Uhat-A||_infinity <= 6 H^2 u. + +Let Ahat=Lhat Uhat. Its row dominance is at least 1/2-6H^2u>1/4, so ||Ahat^-1||_infinity<=4. Also ||Uhat||_infinity<=2 and ||Lhat||_infinity<=1+8h<=9H. Hence + + ||Lhat^-1||_infinity <= ||Uhat|| ||Ahat^-1|| <=8, + ||Uhat^-1||_infinity <= ||Ahat^-1|| ||Lhat|| <=36H. + +The same bounds hold for every computed leading principal pair. They control all forward-substitution prefixes, so no exponential intermediate-growth assumption is being hidden. + +## Rounded solves and total error + +For each row, accumulate its fixed-point dot product exactly with a bounded-word double-precision accumulator and round only the solved coordinate. Forward substitution then produces yhat with + + Lhat yhat = b+e_f, ||e_f||<=u. + +Thus ||yhat||<=8(||b||+u). Back substitution produces xhat with + + Uhat xhat = yhat+e_b, ||e_b||<=2u, + +because rounding a solved coordinate by at most u contributes at most |uhat_ii|u<=2u to that equation. Hence + + Ahat xhat = b+e_f+Lhat e_b, + ||e_f+Lhat e_b|| <=19H u, + ||xhat|| <=4(||b||+19H u). + +Compare with the true solution x=A^-1 b: + + ||xhat-x|| <=2[6H^2u*4(||b||+19H u)+19H u]. + +For u<=1/(1000H^2), this is bounded by the stated conservative 1000H^4u(||b||+1). All partial dot products have at most H terms, factors bounded by 8, y bounded as above and x bounded by the same argument on upper triangular suffixes or their solved equations. An additional O(log H+log(||b||+1)+1) integral bits therefore covers every exact sum and product; the exact rational complex division uses a pivot modulus at least 1/4. Fixed extra limbs hold products, accumulators and divisions. Every scalar update consequently costs O(1) bounded-word instructions. + +## Gaussian caller and credit + +If the actual Gaussian inverse map N has a rigorously truncated bandwidth h=O(sqrt(p)/alpha), truncation norm error eta and ||N-I||<=rho<1/2, apply this theorem to the stored truncated matrix A. Perturbing its right-hand solve changes the ideal inverse by at most ||N^-1||*eta*||A^-1|| times the right-hand norm, so tail and coefficient errors can be invoiced explicitly. Matrix formation, proof of the bandwidth and row-sum tail bound remain the caller's responsibility. This theorem does not assert those facts merely from the name 'Gaussian'. + +The elimination and diagonal-dominance arguments are classical numerical linear algebra. This document supplies a direct native sufficient algorithm and explicit word/precision/cost accounting; it does not claim a new best multiplication exponent or worldwide novelty. The Gaussian resampling caller is credited to the upstream construction and its cited Harvey-Hoeven sources. + + +## The actual periodic Gaussian inverse, with a constructive reduction + +Here the band hypothesis can be proved from the actual upstream formula, rather than postulated. Let s=2, r=t/s>1, q_i=floor(ri+1/2) and beta_i=ri-q_i in [-1/2,1/2]. The source defines N=C T D. Its entries are + + N_ij = sum_(ell congruent j mod s) + exp[-pi alpha^2 (r(ell-i)+beta_i)^2 + pi alpha^2 beta_j^2]. + +For d=|ell-i|>=1, each positive summand is at most + + exp[-pi alpha^2(d^2-d)], + +because |r(ell-i)+beta_i|>=d-1/2 and |beta_j|<=1/2. Truncate the sum to |ell-i|<=h. The omitted absolute row sum is at most + + 2 sum_(d>=h+1) exp[-pi alpha^2 d(d-1)] + <=4 exp[-pi alpha^2 h(h+1)]. + +The final inequality follows by bounding the successive-term ratio by exp[-2pi alpha^2(h+1)]<1/2. Choosing + + h = ceil(sqrt((P ln 2+ln 16)/(pi alpha^2)))+1 + +therefore makes the row-sum tail less than 2^(-P)/4. This is the required h=O(1+sqrt(P)/alpha), with no unknown localization assumption. For s>2h, each retained entry has at most one integer representative in the window, so formation takes O(s h) explicitly generated coefficients. In the actual assembly s is superpolynomial in the precision, hence this condition eventually holds. Smaller s may instead use dense O(s^3) elimination; when s<=2h this is O(s h^2). + +Suppose the actual source's proved estimate ||N-I||<=2^(-alpha^2 theta) is at most 1/8, as it eventually is in the intended regime. All entries are nonnegative real numbers and diagonal retained entries include 1. Removing the positive tail can only decrease the row norm of N-I. Form the retained periodic matrix A to coefficient row-sum error at most 1/32, by allocating O(log h) extra fractional bits per entry. Then ||A-I||<=1/4 and ||A^-1||<=4/3. Its difference from N is the explicitly controlled tail and coefficient-formation error; allocate both to the requested P-bit accuracy instead of fixing the error at 1/32 when final accuracy is required. + +To handle the wraparound without assuming it is an ordinary band matrix, let A0 be A with only its wraparound corner blocks removed. Those off-diagonal entries are nonnegative. With coefficient errors included conservatively, ||A0-I||<=1/4. The removed corner matrix can be written U V, where V selects the first and last h coordinates, and U has those 2h columns equal to the removed corner columns. Thus rank(UV)<=2h, ||V||<=1 and ||U||<=1/4. All these are literal stored matrices. + +Factor the ordinary band matrix A0 using the theorem above. Solve its system for b and for the 2h columns of U. Denote the results by y=A0^-1 b and X=A0^-1 U. The small matrix + + S=I+V X + +has ||S-I||<=1/3, because ||A0^-1||<=4/3. It is itself strictly row diagonally dominant, with inverse norm at most 3/2. Factor this explicit 2h by 2h matrix with the same constructive elimination, solve S z=V y, and return x=y-Xz. Direct multiplication verifies A x=b. Thus the algorithm is not hiding a small inverse oracle. + +The work is O(s h^2): ordinary band factorization costs O(s h^2), its 2h right-hand solves also cost O(s h^2), and the final dense factorization costs O(h^3)<=O(s h^2). Formation of Xz and every small-matrix operation is included. The storage is O(s h). All exact inverse/operator norms above are fixed constants; numerical approximations to X can use per-entry error budget epsilon/(2h), and all finite sums use O(log s+log h) extra guard bits. The explicit band error theorem then gives any requested final accuracy with P+O(log s+log h+log(||b||+1)) bits. This is O(log n) in the intended regime. Bounded-word coefficient preparation is provided by NATIVE-EXP-TABLE.md, including the charged setup. + +Finally, ||N^-1|| is bounded by 8/7 from the original 1/8 estimate. The identity N^-1-A^-1=N^-1(A-N)A^-1 controls truncation and formation errors directly by at most (8/7)(4/3)||A-N||. This completes a positive direct native algorithm for the actual Gaussian inverse, not merely a conditional band-matrix interface. The exact source localization, contraction and resampling identities remain predecessor facts; the native band/Woodbury evaluation and its accounting are the assembly given here. diff --git a/research/pr4554-banded-inverse-check/README.md b/research/pr4554-banded-inverse-check/README.md new file mode 100644 index 0000000000..896b1274ad --- /dev/null +++ b/research/pr4554-banded-inverse-check/README.md @@ -0,0 +1,57 @@ +# Independent check of #4554: native banded inverse and periodic Gaussian Woodbury reduction + +**Task:** independent-moment (independent source-bound coverage). **Base:** canonical main `f626db5c`. **No new κ.** + +Dugongue's #4554 (head `4f46c6a9`) replaces the Gaussian adapter arithmetic with native Word RAM constructions. This package checks one component, `NATIVE-BANDED-INVERSE.md`, against its text, vendored as `PR4554-NATIVE-BANDED-INVERSE.md` (blob `8be8f992`, SHA256-pinned). The component has three parts: +- a pivot-free band elimination with an explicit rounding error bound `1000H⁴u(‖b‖+1)`; +- a proved Gaussian tail and bandwidth choice; +- a rank-2h Woodbury reduction for the periodic wraparound. + +Nothing is imported from #4554's `band_solve.py`, `periodic_solve.py` or `gaussian_inverse.py`. + +## Results + +All checked claims are confirmed. + +**A. Band theorem.** Matrix arithmetic is exact over the Gaussian rationals; norms use complex absolute values. + +The test matrices are near-identity complex band matrices with `ρ = 1/2`: 15 random matrices plus 5 adversarial real ones with every entry of `E` negative. Shapes are `(s, h) ∈ {(12,0), (16,1), (20,2), (24,3), (18,4)}`. + +- **Exact elimination without pivoting.** Pivots are ≥ 1/2, surviving row sums ≤ 3/2, and multipliers ≤ 3. The solve satisfies `‖x‖ ≤ 2‖b‖`. +- **Rounded elimination: 60 solves.** Every multiplier, Schur entry and solved coordinate is rounded once with complex error `< u`, and dot products are accumulated exactly. Three values of `u` are used: the limit `1/(1000H²)`, `2⁻¹⁴` and `2⁻³⁰`. In every solve: + - row sums stay `< 2`, pivots `≥ 1/4` and multipliers `≤ 8`; + - `‖L̂Û − A‖ ≤ 6H²u`, `‖L̂‖ ≤ 9H` and `‖Û‖ ≤ 2`; + - `‖ŷ‖ ≤ 8(‖b‖+u)` and `‖x̂‖ ≤ 4(‖b‖+19Hu)`; + - the final error is at most `1000H⁴u(‖b‖+1)`. The worst observed ratio of error to bound is 8.8·10⁻⁵, consistent with the stated deliberately loose constants. + +**B. Periodic Gaussian inverse.** +- **Summand and tail.** The per-summand distance `|r(ℓ−i)+β_i| ≥ d − 1/2` holds for `r > 1` (2,000 random cases). The tail `2Σ_{d>h} exp(−πα²d(d−1)) ≤ 4exp(−πα²h(h+1))` holds, in log domain, for `α ∈ [2, 6]` and `h ≤ 8`. +- **Choice of `h`.** The stated `h` gives a tail `< 2^(−P)/4` for `P ≤ 4096`. Together these make 5,346 checks. +- **Actual instances.** On 7 instances of `N` with `‖N − I‖ ≤ 1/8` (as low as `r = 1.1875`, where `‖N−I‖ ≈ 9·10⁻³`), using `h = 2`: + - the row tails meet the bound; + - `‖A − I‖ ≤ 1/4`; + - every wraparound entry lies in the first or last h columns, so `A = A0 + UV`; + - `‖A0 − I‖ ≤ 1/4`, `‖U‖ ≤ 1/4` and `‖S − I‖ ≤ 1/3`; + - the Woodbury solve `x = y − Xz` satisfies `A x = b` **exactly** in rational arithmetic. + +**C. Binding and replay.** 8 statements are bound verbatim to the vendored text. I also replayed #4554's registered package at its head: `PASS_REGISTERED_RESEARCH_PACKAGE research/native-word-fft outcome=completed` in 16.6 s. + +## Verify + +```text +python3 -B research/pr4554-banded-inverse-check/verify.py +``` + +This runs in about 16 s with the standard library only. + +## Limitations + +- **Out of scope:** the rest of #4554 (roots and recovery, phase and exp tables, native stability, outer assembly), and the source estimate `‖N − I‖ ≤ 2^(−α²θ)`, which is a predecessor fact. This package only checks that its test instances lie in that regime. +- **Floating-point norms.** Instance entries of `N` and all norm comparisons use floats. Matrix arithmetic and the Woodbury identity are exact. +- **Finite checks.** The universal statements rest on #4554's written induction, which these finite checks corroborate. + +## Open next questions + +- An independent check of `NATIVE-STABILITY.md` or `ROOTS-RECOVERY.md`, the other precision-carrying components of #4554. First step: an exact-rational replay of their stated per-step error recurrences. + +Credit: Dugongue (#4554); the Gaussian resampling construction and the Harvey–van der Hoeven sources it cites. Prepared by Ryanson Jonathan with substantial Claude Code (Claude Opus 5.5) assistance. diff --git a/research/pr4554-banded-inverse-check/SOURCES.json b/research/pr4554-banded-inverse-check/SOURCES.json new file mode 100644 index 0000000000..a2a5ea8bbd --- /dev/null +++ b/research/pr4554-banded-inverse-check/SOURCES.json @@ -0,0 +1,18 @@ +{ + "schema_version": 1, + "base_commit": "f626db5c3a3d781d8f7695038776c621651f75e1", + "files": [ + { + "path": "certificates/selected-result.json", + "usage": "provenance", + "git_blob": "ed33caef5c5827693e9572d252a66f99cbb73d40", + "sha256": "5732707cc1e460b090b53f2463162749ad860de878a18f7408938671d295a9d9" + }, + { + "path": "docs/research/agent-task-catalog.json", + "usage": "provenance", + "git_blob": "64bfc5af6709ed9811171e0a2108f05e1219cb65", + "sha256": "352a6184713b8843d96fd2c2071eb7aeaa655906b7c51d45ad444912666c8a96" + } + ] +} diff --git a/research/pr4554-banded-inverse-check/TASK.json b/research/pr4554-banded-inverse-check/TASK.json new file mode 100644 index 0000000000..c4e5f721f9 --- /dev/null +++ b/research/pr4554-banded-inverse-check/TASK.json @@ -0,0 +1,100 @@ +{ + "schema_version": 1, + "task_id": "independent-moment", + "priority_revision": 7, + "base_commit": "f626db5c3a3d781d8f7695038776c621651f75e1", + "original_card": { + "id": "independent-moment", + "question": "Does a new paid-profile or mixture-ceiling claim have independent source-bound coverage of its actual changed ledger?", + "starting_points": [ + "docs/research/research-family-status.md", + "certificates/selected-result.json", + "research/staged544-mixed-fourth-entry/README.md", + "research/staged-fourth-mixture-ceiling/README.md", + "research/r437-two-window-carrier/README.md" + ], + "methods": [ + "Independently reconstruct one changed paid ledger or pooled dual with exact source and model binding.", + "Check a genuinely uncovered moment implication with a different enclosure method and semantic controls." + ], + "acceptance_gate": "Demonstrate independence, source binding, positive cases and semantic failure controls. Name the existing checks and the missing coverage. More digits or another equivalent evaluation of an unchanged root is not by itself new verification coverage.", + "first_step": { + "command": [ + "python3", + "-B", + "scripts/research_baseline.py" + ], + "expected": "PASS_RESEARCH_BASELINE_BINDINGS_ONLY", + "inspect": "Choose one uncovered implication in #4088 literal mixed fill, #4091 pooled dual/model binding or #4067 changed frame/host bill. #4091 pins an earlier #4088 head: validate word, capacity and stock equivalence before applying that ceiling to the later fill.", + "scope": "Source-binding orientation only; no scientific replay, theorem acceptance or supplier admission. Budget any heavy replay separately." + }, + "first_experiment": "Independently reconstruct one changed ledger or pooled certificate with a positive case and a semantic failure control. Compare exact sources and name what existing checks miss. Keep this separate from supplier-boundary-audit (ordinary/outer dispatch) and uniform-transfer-audit (all-size proof). Stop after the chosen implication; more digits of an unchanged root are not additional mechanism coverage.", + "bottleneck_link": "Independent pricing prevents unsupported frontier claims and can reveal genuine accounting slack.", + "allowed_variations": [ + "Reconstruct either a source ledger or a numeric enclosure with a different method." + ], + "scope_change_triggers": [ + "Importing the primary function and calling its output independent verification.", + "Calling an aggregate histogram or source archive an independently accepted physical supplier." + ], + "overlap_terms": [ + "independent moment", + "rank histogram", + "exact enclosure", + "2979", + "rank7 joint endframes", + "paired extraction", + "joint lowering", + "4088", + "4091", + "4067", + "staged544-mixed-fourth-entry", + "staged-fourth-mixture-ceiling", + "r437-two-window-carrier" + ], + "existing_checks": [ + { + "scope": "PR #4088: mixed local/owner/phase ledger and moment reported.", + "evidence": "research/staged544-mixed-fourth-entry/README.md", + "remaining": "Independent binding of changed construction to paid histogram, distinct from full admission.", + "status": "author-checked", + "source_commit": "72d7efc181edb21806a878f174847c96dd410ab0", + "selected_record_sha256": "bb2eaa58647b353f88d0aec2eafb9d30df1be53589231b23e12ceef3e4bad129" + }, + { + "scope": "PR #4091: sharp pooled ceiling proposed for fixed THREE/FOUR words and inventory.", + "evidence": "research/staged-fourth-mixture-ceiling/README.md", + "remaining": "Validate dual, exhaustive pattern scope and exact source equivalence before applying to newer fills.", + "status": "author-checked", + "source_commit": "a29bdd8f56cd18b805917775cdec09ff935eddc4", + "selected_record_sha256": "bb2eaa58647b353f88d0aec2eafb9d30df1be53589231b23e12ceef3e4bad129" + }, + { + "scope": "PR #4067: paid (2,2) to (1,3) exchange and fixed OQ4 host reported.", + "evidence": "research/r437-two-window-carrier/README.md", + "remaining": "Independent full cost delta; remaining caller/precision/stock obligations are separate.", + "status": "author-checked", + "source_commit": "a95f11c52e3ca1e4c5171b3c20d8f4b0815aee2e", + "selected_record_sha256": "bb2eaa58647b353f88d0aec2eafb9d30df1be53589231b23e12ceef3e4bad129" + } + ] + }, + "scope_status": "unchanged", + "actual_scope": "Does a new paid-profile or mixture-ceiling claim have independent source-bound coverage of its actual changed ledger?", + "scope_change_reason": null, + "bottleneck_link": "Independent coverage of the native Gaussian-inverse adapter proposed for the Word RAM port of the multiplication route.", + "budget_minutes": 45, + "stop_rule": "Stop at NATIVE-BANDED-INVERSE.md.", + "outcome": "completed", + "evidence": "Exact Gaussian-rational elimination on 20 near-identity band matrices (incl. adversarial rho=1/2): exact and rounded invariants and the 1000H^4u(||b||+1) bound in 60 rounded solves (worst ratio 8.8e-5); 5,346 tail/h checks; 7 actual periodic Gaussian instances with exact Woodbury solves A x = b; 8 statements bound; #4554 package replays.", + "author": "Ryanson Jonathan, with Claude Code (Claude Opus 5.5) assistance", + "ai_assistance": "Substantial Claude Code (Claude Opus 5.5) assistance; no independent human review.", + "open_questions": [ + { + "question": "Independent check of #4554 NATIVE-STABILITY.md or ROOTS-RECOVERY.md.", + "why_it_matters": "Other precision-carrying components of the native adapter.", + "first_step": "Exact-rational replay of their per-step error recurrences.", + "blocked_by": "None." + } + ] +} diff --git a/research/pr4554-banded-inverse-check/VERIFY.json b/research/pr4554-banded-inverse-check/VERIFY.json new file mode 100644 index 0000000000..673b0cfc74 --- /dev/null +++ b/research/pr4554-banded-inverse-check/VERIFY.json @@ -0,0 +1,18 @@ +{ + "schema_version": 1, + "status": "ready", + "python": "3.12", + "requirements": [], + "inputs": [], + "commands": [ + { + "argv": [ + "{python}", + "-B", + "{package}/verify.py" + ], + "timeout_seconds": 300, + "checks": "Exact band elimination invariants and rounded error bound; Gaussian tail and h choice; periodic corner split and exact Woodbury solves on actual N instances; binding to vendored spec (sha256)." + } + ] +} diff --git a/research/pr4554-banded-inverse-check/verify.py b/research/pr4554-banded-inverse-check/verify.py new file mode 100644 index 0000000000..6b39034f2a --- /dev/null +++ b/research/pr4554-banded-inverse-check/verify.py @@ -0,0 +1,316 @@ +#!/usr/bin/env python3 +"""Independent check of PR #4554 (Dugongue) NATIVE-BANDED-INVERSE.md, written from the vendored text only. + +Nothing is imported from #4554's band_solve.py, periodic_solve.py or gaussian_inverse.py. + A. Band theorem. Exact Gaussian-rational elimination without pivoting on random and adversarial (rho = 1/2) + complex near-identity band matrices: pivots >= 1/2, row sums <= 3/2, multipliers <= 3, ||A^-1|| <= 2. + Rounded elimination with explicit rounding (complex error < u, every multiplier, Schur entry and solved + coordinate rounded once; dot products accumulated exactly): row sums < 2, pivots >= 1/4, multipliers <= 8, + ||Lhat Uhat - A|| <= 6H^2 u, ||Lhat|| <= 9H, ||Uhat|| <= 2, ||yhat|| <= 8(||b||+u), ||xhat|| <= 4(||b||+19Hu), + and the final ||xhat - x|| <= 1000 H^4 u (||b||+1), for u from 1/(1000H^2) down to 2^-30. + B. Periodic Gaussian inverse. The per-summand bound |r(l-i)+beta_i| >= d - 1/2; the tail + 2 sum_{d>h} exp(-pi a^2 d(d-1)) <= 4 exp(-pi a^2 h(h+1)) (log domain); the choice of h against 2^-P/4; then on + actual instances of N (s > 2h, ||N - I|| <= 1/8): row tails, ||A - I|| <= 1/4, the corner split A = A0 + UV + with ||A0 - I|| <= 1/4, ||U|| <= 1/4, ||S - I|| <= 1/3, and the Woodbury solve x = y - Xz satisfying A x = b + exactly in rational arithmetic. +Norms are infinity norms with complex absolute values (floats, for comparisons with stated slack only); all matrix +arithmetic is exact. Ryanson Jonathan with Claude Code (Claude Opus 5.5) assistance. Apache-2.0. +""" +from fractions import Fraction as Fr +import hashlib +import json +import math +from pathlib import Path +import random +import sys + +HERE = Path(__file__).resolve().parent +SPEC = HERE / "PR4554-NATIVE-BANDED-INVERSE.md" # #4554 head 4f46c6a9, research/native-word-fft/NATIVE-BANDED-INVERSE.md, blob 8be8f992 +SPEC_SHA256 = "1259b4242ebf6ac7f14bfe72e371afb8dc530b2755d1fd82228d9dc6f01634ba" +ZERO = (Fr(0), Fr(0)) +ONE = (Fr(1), Fr(0)) + + +def need(condition, message): + if not condition: + raise AssertionError(message) + + +def cadd(a, b): + return (a[0] + b[0], a[1] + b[1]) + + +def csub(a, b): + return (a[0] - b[0], a[1] - b[1]) + + +def cmul(a, b): + return (a[0] * b[0] - a[1] * b[1], a[0] * b[1] + a[1] * b[0]) + + +def cdiv(a, b): + den = b[0] * b[0] + b[1] * b[1] + return ((a[0] * b[0] + a[1] * b[1]) / den, (a[1] * b[0] - a[0] * b[1]) / den) + + +def cabs(a): + return math.hypot(float(a[0]), float(a[1])) + + +def rnd(a, u): + """Round each part to the grid u/2: complex error at most u*sqrt(2)/4 < u.""" + g = u / 2 + return (Fr(round(a[0] / g)) * g, Fr(round(a[1] / g)) * g) + + +def norm_mat(M): + return max(sum(cabs(x) for x in row) for row in M) + + +def norm_vec(v): + return max(cabs(x) for x in v) + + +def matvec(M, x): + out = [] + for row in M: + acc = ZERO + for a, b in zip(row, x): + if a != ZERO: + acc = cadd(acc, cmul(a, b)) + out.append(acc) + return out + + +def random_band(rng, s, h, rho, adversarial): + A = [[ZERO] * s for _ in range(s)] + for i in range(s): + cols = [j for j in range(max(0, i - h), min(s, i + h + 1))] + raw = [] + for j in cols: + if adversarial: + z = (-1.0, 0.0) if j != i else (-1.0, 0.0) + else: + ang = rng.uniform(0, 2 * math.pi) + z = (math.cos(ang) * rng.random(), math.sin(ang) * rng.random()) + raw.append(z) + tot = sum(math.hypot(*z) for z in raw) + scale = rho * 0.999999 / tot if tot else 0 + for j, z in zip(cols, raw): + e = (Fr(z[0] * scale).limit_denominator(1 << 30), Fr(z[1] * scale).limit_denominator(1 << 30)) + A[i][j] = cadd(ONE, e) if i == j else e + need(sum(cabs(A[i][j] if j != i else csub(A[i][j], ONE)) for j in cols) <= rho, "||E|| <= rho") + return A + + +def exact_elimination(A, h): + s = len(A) + M = [row[:] for row in A] + worst = {"pivot_min": 9.0, "rowsum_max": 0.0, "mult_max": 0.0} + for k in range(s): + piv = M[k][k] + worst["pivot_min"] = min(worst["pivot_min"], cabs(piv)) + for i in range(k + 1, min(s, k + h + 1)): + if M[i][k] == ZERO: + continue + l = cdiv(M[i][k], piv) + worst["mult_max"] = max(worst["mult_max"], cabs(l)) + for j in range(k + 1, min(s, k + h + 1)): + M[i][j] = csub(M[i][j], cmul(l, M[k][j])) + M[i][k] = ZERO + worst["rowsum_max"] = max(worst["rowsum_max"], sum(cabs(x) for x in M[i][k + 1:])) + return worst + + +def exact_solve(A, b): + s = len(A) + M = [row[:] + [b[i]] for i, row in enumerate(A)] + for k in range(s): + for i in range(k + 1, s): + if M[i][k] != ZERO: + l = cdiv(M[i][k], M[k][k]) + M[i] = [csub(x, cmul(l, y)) for x, y in zip(M[i], M[k])] + x = [ZERO] * s + for i in range(s - 1, -1, -1): + acc = M[i][s] + for j in range(i + 1, s): + acc = csub(acc, cmul(M[i][j], x[j])) + x[i] = cdiv(acc, M[i][i]) + return x + + +def rounded_solve(A, b, h, u): + s = len(A) + H = h + 1 + U = [row[:] for row in A] + L = [[ZERO] * s for _ in range(s)] + st = {"pivot_min": 9.0, "mult_max": 0.0, "rowsum_max": 0.0} + for k in range(s): + piv = U[k][k] + st["pivot_min"] = min(st["pivot_min"], cabs(piv)) + for i in range(k + 1, min(s, k + h + 1)): + l = rnd(cdiv(U[i][k], piv), u) + st["mult_max"] = max(st["mult_max"], cabs(l)) + L[i][k] = l + for j in range(k + 1, min(s, k + h + 1)): + U[i][j] = rnd(csub(U[i][j], cmul(l, U[k][j])), u) + U[i][k] = ZERO + st["rowsum_max"] = max(st["rowsum_max"], sum(cabs(x) for x in U[i][k + 1:])) + for i in range(s): + L[i][i] = ONE + LU = [[ZERO] * s for _ in range(s)] + for i in range(s): + for j in range(s): + acc = ZERO + for k in range(max(0, i - h, j - h), min(i, j) + 1): + acc = cadd(acc, cmul(L[i][k], U[k][j])) + LU[i][j] = acc + resid = norm_mat([[csub(LU[i][j], A[i][j]) for j in range(s)] for i in range(s)]) + y = [ZERO] * s + for i in range(s): + acc = b[i] + for k in range(max(0, i - h), i): + acc = csub(acc, cmul(L[i][k], y[k])) + y[i] = rnd(acc, u) + x = [ZERO] * s + for i in range(s - 1, -1, -1): + acc = y[i] + for j in range(i + 1, min(s, i + h + 1)): + acc = csub(acc, cmul(U[i][j], x[j])) + x[i] = rnd(cdiv(acc, U[i][i]), u) + st.update(resid=resid, Lnorm=norm_mat(L), Unorm=norm_mat(U), ynorm=norm_vec(y), xnorm=norm_vec(x)) + return x, st + + +def check_band(rng): + runs, worst_ratio = 0, 0.0 + for s, h in ((12, 0), (16, 1), (20, 2), (24, 3), (18, 4)): + H = h + 1 + for adversarial in (False, True): + for trial in range(3 if not adversarial else 1): + A = random_band(rng, s, h, Fr(1, 2), adversarial) + ex = exact_elimination(A, h) + need(ex["pivot_min"] >= 0.5 - 1e-12 and ex["rowsum_max"] <= 1.5 + 1e-12 and ex["mult_max"] <= 3 + 1e-12, + "exact pivots, row sums, multipliers") + b = [(Fr(rng.randint(-9, 9), 4), Fr(rng.randint(-9, 9), 4)) for _ in range(s)] + x = exact_solve(A, b) + need(norm_vec(x) <= 2 * norm_vec(b) + 1e-12, "||A^-1|| <= 2 on b") + for u in (Fr(1, 1000 * H * H), Fr(1, 2 ** 14), Fr(1, 2 ** 30)): + xh, st = rounded_solve(A, b, h, u) + uf, nb = float(u), norm_vec(b) + need(st["rowsum_max"] < 2 and st["pivot_min"] >= 0.25 and st["mult_max"] <= 8, "rounded invariants") + need(st["resid"] <= 6 * H * H * uf and st["Lnorm"] <= 9 * H and st["Unorm"] <= 2, "factor bounds") + need(st["ynorm"] <= 8 * (nb + uf) and st["xnorm"] <= 4 * (nb + 19 * H * uf), "solve growth") + err = norm_vec([csub(p, q) for p, q in zip(xh, x)]) + bound = 1000 * H ** 4 * uf * (nb + 1) + need(err <= bound, "final error bound") + worst_ratio = max(worst_ratio, err / bound) + runs += 1 + return runs, worst_ratio + + +def check_gaussian_tail(): + rng = random.Random(7) + for _ in range(2000): # |r(l-i) + beta_i| >= d - 1/2 for r > 1, |beta| <= 1/2 + r = 1 + rng.random() * 3 + d = rng.randint(1, 20) * rng.choice((-1, 1)) + beta = rng.uniform(-0.5, 0.5) + need(abs(r * d + beta) >= abs(d) - 0.5, "per-summand distance") + tails = 0 + for a10 in range(20, 61, 5): + a = a10 / 10 + for h in range(1, 9): + logs = [-math.pi * a * a * d * (d - 1) for d in range(h + 1, h + 60)] + m = max(logs) + lhs = math.log(2) + m + math.log(sum(math.exp(v - m) for v in logs)) + rhs = math.log(4) - math.pi * a * a * h * (h + 1) + need(lhs <= rhs and math.exp(-2 * math.pi * a * a * (h + 1)) < 0.5, "tail bound") + tails += 1 + for P in range(1, 4097, 7): + h = math.ceil(math.sqrt((P * math.log(2) + math.log(16)) / (math.pi * a * a))) + 1 + need(math.pi * a * a * h * (h + 1) > P * math.log(2) + math.log(16), "h choice gives tail < 2^-P/4") + tails += 1 + return tails + + +def gaussian_N(s, t, a): + r = t / s + beta = [r * i - math.floor(r * i + 0.5) for i in range(s)] + N = [[0.0] * s for _ in range(s)] + for i in range(s): + for l in range(i - 14, i + 15): + j = l % s + N[i][j] += math.exp(-math.pi * a * a * (r * (l - i) + beta[i]) ** 2 + math.pi * a * a * beta[j] ** 2) + return N, beta + + +def check_periodic(rng): + instances = 0 + for s, t, a in ((24, 37, 2.0), (30, 47, 2.0), (40, 61, 2.5), (33, 70, 2.0), (28, 45, 3.0), (32, 38, 2.0), (40, 50, 2.0)): + N, beta = gaussian_N(s, t, a) + nI = max(sum(abs(N[i][j] - (1 if i == j else 0)) for j in range(s)) for i in range(s)) + need(nI <= 0.125, "instance in the stated regime ||N - I|| <= 1/8") + h = 2 + need(s > 2 * h, "s > 2h") + r = t / s + A = [[ZERO] * s for _ in range(s)] + tail_max = 0.0 + for i in range(s): + full = sum(N[i]) + kept = 0.0 + for l in range(i - h, i + h + 1): + j = l % s + v = math.exp(-math.pi * a * a * (r * (l - i) + beta[i]) ** 2 + math.pi * a * a * beta[j] ** 2) + A[i][j] = (Fr(v), Fr(0)) + kept += v + tail_max = max(tail_max, full - kept) + need(tail_max <= 4 * math.exp(-math.pi * a * a * h * (h + 1)) + 1e-15, "row tail") + need(norm_mat([[csub(A[i][j], ONE if i == j else ZERO) for j in range(s)] for i in range(s)]) <= 0.25, "||A-I||") + # Corner split: A0 keeps |i-j| <= h; the corner part has columns only in the first and last h coordinates. + A0 = [[A[i][j] if abs(i - j) <= h else ZERO for j in range(s)] for i in range(s)] + corner_cols = list(range(h)) + list(range(s - h, s)) + Ucols = [[A[i][j] if abs(i - j) > h else ZERO for j in corner_cols] for i in range(s)] + for i in range(s): + for j in range(s): + if abs(i - j) > h and A[i][j] != ZERO: + need(j in corner_cols, "corner entries lie in the first/last h columns") + need(norm_mat([[csub(A0[i][j], ONE if i == j else ZERO) for j in range(s)] for i in range(s)]) <= 0.25, "||A0-I||") + need(norm_mat(Ucols) <= 0.25, "||U|| <= 1/4") + b = [(Fr(rng.randint(-5, 5)), Fr(rng.randint(-5, 5))) for _ in range(s)] + y = exact_solve(A0, b) + Xcols = [exact_solve(A0, [Ucols[i][c] for i in range(s)]) for c in range(2 * h)] + X = [[Xcols[c][i] for c in range(2 * h)] for i in range(s)] + VX = [X[j] for j in corner_cols] + S = [[cadd(VX[p][q], ONE if p == q else ZERO) for q in range(2 * h)] for p in range(2 * h)] + need(norm_mat([[VX[p][q] for q in range(2 * h)] for p in range(2 * h)]) <= 1 / 3, "||S - I|| <= 1/3") + z = exact_solve(S, [y[j] for j in corner_cols]) + Xz = matvec(X, z) + x = [csub(y[i], Xz[i]) for i in range(s)] + need(matvec(A, x) == b, "Woodbury solve: A x = b exactly") + instances += 1 + return instances + + +def check_binding(): + raw = SPEC.read_bytes() + need(hashlib.sha256(raw).hexdigest() == SPEC_SHA256, "vendored NATIVE-BANDED-INVERSE.md bytes") + text = raw.decode("utf-8") + for s in ("1000 H^4 u (||b||_infinity+1),", "||Lhat Uhat-A||_infinity <= 6 H^2 u.", "||e_f+Lhat e_b|| <=19H u,", + "||xhat|| <=4(||b||+19H u).", "<=4 exp[-pi alpha^2 h(h+1)].", + "h = ceil(sqrt((P ln 2+ln 16)/(pi alpha^2)))+1", "S=I+V X", "return x=y-Xz"): + need(s in text, f"spec states: {s}") + return 8 + + +def main(): + rng = random.Random(4554) + runs, ratio = check_band(rng) + report = {"status": "PASS_PR4554_BANDED_INVERSE_CHECK", "rounded_band_solves": runs, + "worst_error_over_bound": f"{ratio:.3e}", "gaussian_tail_and_h_checks": check_gaussian_tail(), + "periodic_woodbury_instances": check_periodic(rng), "spec_statements_bound": check_binding(), + "new_kappa_claimed": False} + print(json.dumps(report, indent=1)) + + +if __name__ == "__main__": + sys.exit(main())