From ed4a13f4d4b2b55b06d21bfbb4deedcfdaa0aecb Mon Sep 17 00:00:00 2001 From: Mingli Yuan Date: Thu, 27 Aug 2026 18:05:30 +0800 Subject: [PATCH] research: compile AMP polynomial and matrix layers --- sonnet/README.md | 8 + .../00-problem-frontier.md | 112 +++++ .../01-sparse-compiler-theorems.md | 180 ++++++++ .../02-benchmark-results.md | 118 +++++ .../03-disposition.md | 87 ++++ .../amp-polynomial-matrix-compiler/README.md | 30 ++ .../amp_escape_compiler.py | 417 ++++++++++++++++++ .../test_amp_polynomial_matrix_compiler.py | 180 ++++++++ 8 files changed, 1132 insertions(+) create mode 100644 sonnet/amp-polynomial-matrix-compiler/00-problem-frontier.md create mode 100644 sonnet/amp-polynomial-matrix-compiler/01-sparse-compiler-theorems.md create mode 100644 sonnet/amp-polynomial-matrix-compiler/02-benchmark-results.md create mode 100644 sonnet/amp-polynomial-matrix-compiler/03-disposition.md create mode 100644 sonnet/amp-polynomial-matrix-compiler/README.md create mode 100644 sonnet/amp-polynomial-matrix-compiler/amp_escape_compiler.py create mode 100644 tests/research/test_amp_polynomial_matrix_compiler.py diff --git a/sonnet/README.md b/sonnet/README.md index f3614749..8f0f3fc4 100644 --- a/sonnet/README.md +++ b/sonnet/README.md @@ -190,6 +190,14 @@ tail outside the finite carrier. The AMP line earns `EXPAND`; every frozen workload eliminates a surreal runtime, so the overall gate remains `NARROW` until an interacting residual yields a measured computation advantage. +[`amp-polynomial-matrix-compiler/`](amp-polynomial-matrix-compiler/) performs +that first interacting-residual gate on `x -> x^d+t`. The AMP exponential-ray +basis turns composition into an exact sparse nilpotent matrix and the +long-horizon coordinate into a triangular linear solve. It eliminates +exponential symbolic support and supplies exact residual certificates, while a +strong logarithmic recurrence prevents any claim of universal numerical +speedup. The result is `EXPAND-NARROW` and remains Sonnet-local. + ## Research-local calibration — the \(S^6\) complex structure claim [`s6-complex-arithmetic-tower/`](s6-complex-arithmetic-tower/) studies a diff --git a/sonnet/amp-polynomial-matrix-compiler/00-problem-frontier.md b/sonnet/amp-polynomial-matrix-compiler/00-problem-frontier.md new file mode 100644 index 00000000..fbf8f931 --- /dev/null +++ b/sonnet/amp-polynomial-matrix-compiler/00-problem-frontier.md @@ -0,0 +1,112 @@ +# Problem frontier: do AMP polynomials and matrices simplify an algorithm? + +Status: frozen contract for issue +[#152](https://github.com/mountain/process-geometry/issues/152). + +## 1. The question + +Issue #150 identified the finite AMP polynomial-like family + +\[ +\sum_{\gamma,n}^{\mathrm{finite}} +c_{\gamma,n}x^\gamma(\log x)^n +\] + +and the sparse operator matrices induced by `A`, `M`, and `P` on the degree +lattice `(gamma,n)`. Algebraic closure alone does not establish algorithmic +value. This phase asks whether the representation reduces a frozen task under +same-information baselines. + +## 2. Frozen power-dominant task + +Take + +\[ +f(x)=x^d+t, +\qquad d\ge2, +\qquad t>0, +\] + +on a declared positive chart near infinity. Put + +\[ +y=\log x, +\qquad +q=e^{-y}. +\] + +Then + +\[ +F(y)=dy+\log(1+tq^d), +\qquad +g(q)=e^{-F(y)}=\frac{q^d}{1+tq^d}. +\] + +The main observer is the normalized long-horizon quantity + +\[ +G_N(y)=d^{-N}F^{\circ N}(y) +\] + +and its limit when the asymptotic coordinate converges. This is deliberately +narrower than reconstructing the full iterate or orbit. + +## 3. Same-information baselines + +Three paths receive separate ledgers. + +1. **Expanded symbolic baseline:** form `f^[N](x)` as one ordinary expanded + polynomial. +2. **Strong numerical baseline:** update `F(y)` directly in the logarithmic + chart, accumulate the normalized correction, and stop when floating-point + correction is zero. It is forbidden to expand the polynomial. +3. **AMP compiler:** compile a finite polynomial-like escape coordinate once + from a sparse substitution matrix, then evaluate it for the declared + observer. + +The expanded baseline measures symbolic support only. It may not be used as +the sole numerical competitor. + +## 4. Metrics + +- exact expanded support count; +- observer order and nonzero polynomial-like terms; +- dense versus sparse matrix entries; +- exact compilation and residual certificate; +- strong-baseline executed steps; +- compile-once/evaluate-many online work; +- numerical error across observer orders; +- chart failure outside the asymptotic domain; +- decoder and output scope. + +## 5. Acceptance and kill conditions + +The representation earns algorithmic credit only if: + +- the basis makes the declared transport exactly sparse; +- the matrix construction changes the solve, not only its notation; +- the result is replayable without hidden symbolic expansion; +- a strong recurrence baseline is reported; +- compilation, online work, and output precision are separated. + +Narrow or stop a claim if: + +- an ordinary Taylor or monomial basis is relabelled AMP without a support + advantage; +- a dense generic eigensolver replaces an available triangular solve; +- symbolic expansion is treated as the only numerical baseline; +- the finite truncation is evaluated outside its chart without a refusal; +- one scalar long-time observer is presented as full-orbit reconstruction; +- a classical Böttcher/Koopman result is claimed as new. + +## 6. Claim ceiling + +This phase does not claim a new Böttcher theorem, generic Koopman solver, +complexity-class improvement, Ising solver, or Public API. + +```text +Epistemic maturity: T1 exact finite compiler + bounded numerical calibration +Engineering status: Sonnet-local Python +Mathematical Core: unchanged +``` diff --git a/sonnet/amp-polynomial-matrix-compiler/01-sparse-compiler-theorems.md b/sonnet/amp-polynomial-matrix-compiler/01-sparse-compiler-theorems.md new file mode 100644 index 00000000..1b1411d8 --- /dev/null +++ b/sonnet/amp-polynomial-matrix-compiler/01-sparse-compiler-theorems.md @@ -0,0 +1,180 @@ +# Sparse AMP compiler theorems + +Status: exact over rational `t>0` at every fixed observer order. + +## 1. Polynomial-like observer basis + +On the `y=log x` chart, use + +\[ +q^k=e^{-ky}=x^{-k}, +\qquad 1\le k\le K. +\] + +The compiled coordinate has the finite AMP form + +\[ +H_K(y)=y+\sum_{k=1}^K h_kq^k. +\] + +It is a finite slice of the exponential-logarithmic algebra: the affine `y` +term records the Power direction, and the `q` ray records the completed +Addition residual around the power-dominant chart. + +## 2. Matrix-like composition + +Let `C_K` be composition by + +\[ +g(q)=\frac{q^d}{1+tq^d} +\] + +on the truncated ray. Its entries are + +\[ +(C_K)_{r,k}=[q^r]g(q)^k. +\] + +Using the negative-binomial expansion, + +\[ +g(q)^k +=q^{dk}(1+tq^d)^{-k} +=\sum_{j\ge0} +(-1)^j\binom{k+j-1}{j}t^j q^{d(k+j)}. +\] + +Therefore + +\[ +(C_K)_{d(k+j),k} +=(-1)^j\binom{k+j-1}{j}t^j, +\] + +and every other entry is zero. + +**Proposition 2.1.** `C_K` is nilpotent. In particular, + +\[ +C_K^r=0 +\quad\text{whenever}\quad +d^r>K. +\] + +**Proof.** One application sends a monomial of degree `k` to degrees at +least `dk`. After `r` applications, every surviving degree is at least +`d^r k`. No positive degree survives the observer cutoff when `d^r>K`. +QED. + +This nilpotence is the finite-observer version of power-driven scale escape. + +## 3. The linear eigenproblem + +Write + +\[ +u(q)=\log(1+tq^d) +=\sum_{j\ge1}\frac{(-1)^{j+1}}{j}t^jq^{dj}. +\] + +Composition gives + +\[ +H_K(F(y)) +=dy+u(q)+\sum_{k=1}^K h_k g(q)^k. +\] + +Thus the truncated conjugacy condition + +\[ +H_K\circ F=dH_K+O(q^{K+1}) +\] + +is exactly + +\[ +(dI-C_K)h=u_{\le K}. +\] + +**Theorem 3.1.** This equation has a unique rational solution for rational +`t`. + +**Proof.** `C_K` strictly raises degree, so `dI-C_K` is triangular with +nonzero diagonal `d`. Equivalently, nilpotence gives the finite inverse + +\[ +(dI-C_K)^{-1} +=\frac1d\sum_{r\ge0}^{\mathrm{finite}} +\left(\frac{C_K}{d}\right)^r. +\] + +All entries remain rational. QED. + +For `d=2`, `t=1`, and `K=10`, the compiler obtains + +\[ +H_{10}(y) +=y+\frac12q^2-\frac13q^6+\frac58q^8-\frac9{10}q^{10}. +\] + +The finite eigenrelation replays exactly. Its first omitted residual is + +\[ +2q^{12}. +\] + +## 4. Long-horizon observer + +The exact, untruncated coordinate is the logarithm of a Böttcher coordinate +and satisfies + +\[ +H(F(y))=dH(y). +\] + +Consequently, + +\[ +d^{-N}H(F^{\circ N}(y))=H(y). +\] + +When the correction `H(z)-z` vanishes along the escaping orbit, + +\[ +\lim_{N\to\infty}d^{-N}F^{\circ N}(y)=H(y). +\] + +The finite compiler approximates this observer directly, without constructing +the degree-`d^N` iterate. + +## 5. Classical boundary + +This construction meets established mathematics: + +- Böttcher coordinates conjugate a degree-`d` polynomial near infinity to + `z -> z^d`; +- the Green/escape-rate function is the logarithmic long-horizon observer; +- the Koopman operator acts linearly on observables by composition; +- Carleman-style methods represent nonlinear composition in an infinite + function basis. + +The AMP-specific question is narrower: does the arithmetic chart select the +right sparse dictionary and typed completion automatically? The present +example gives one positive calibration, not a novelty claim about these +classical structures. + +Primary references used as boundaries: + +1. B. O. Koopman, + [“Hamiltonian Systems and Transformation in Hilbert Space”](https://www.pnas.org/doi/10.1073/pnas.17.5.315), + 1931. +2. C. Favre and T. Gauthier, + [“The arithmetic of polynomial dynamical pairs”](https://arxiv.org/abs/2004.13801), + including formal Böttcher expansions at infinity. +3. L. DeMarco, K. Lindsey, + [“Convergence properties of the Gronwall area formula for quadratic Julia sets”](https://arxiv.org/abs/1405.1933), + including coefficient-level numerical use of Böttcher maps. +4. M. J. Colbrook, + [structure-preserving finite approximations of Koopman operators](https://arxiv.org/abs/2209.02244), + a distinct data-driven setting that reinforces the need to state the + dictionary, truncation, and convergence target. diff --git a/sonnet/amp-polynomial-matrix-compiler/02-benchmark-results.md b/sonnet/amp-polynomial-matrix-compiler/02-benchmark-results.md new file mode 100644 index 00000000..a617f258 --- /dev/null +++ b/sonnet/amp-polynomial-matrix-compiler/02-benchmark-results.md @@ -0,0 +1,118 @@ +# Benchmark results + +Status: exact support/certificate measurements plus bounded floating-point +calibration. + +## 1. Frozen instance + +Use + +\[ +d=2, +\qquad +t=1, +\qquad +y_0=1.5, +\qquad +N=100. +\] + +The strong numerical recurrence gives + +\[ +G_{100}(y_0)=1.5248559772600594. +\] + +It detects that the interaction correction has underflowed to zero after nine +executed steps in binary64 arithmetic. This early stop is retained as a +positive baseline result. + +## 2. Symbolic support + +For positive `t`, the fully expanded iterate has + +\[ +\#\operatorname{supp}(f^{\circ N})=d^{N-1}+1. +\] + +For the frozen instance this is + +\[ +2^{99}+1 +=633825300114114700748351602689 +\] + +ordinary polynomial terms. The AMP compiler never constructs this +polynomial; its observer state remains `K` ray coefficients plus the affine +`y` term. + +This is a real symbolic and storage simplification. It does not imply the +same factor against direct numerical recurrence. + +## 3. Sparse compilation and accuracy + +| Order `K` | Dense entries | Sparse entries | Nonzero `h_k` | First residual | Absolute error | +|---:|---:|---:|---:|---:|---:| +| 6 | 36 | 6 | 2 | `5/4 q^8` | `3.58e-6` | +| 10 | 100 | 15 | 4 | `2 q^12` | `1.49e-8` | +| 14 | 196 | 28 | 6 | `-3 q^16` | `4.82e-11` | +| 20 | 400 | 55 | 9 | `144/11 q^22` | `3.13e-14` | +| 30 | 900 | 120 | 14 | `-2523/16 q^32` | below binary64 distinction | + +All coefficients, sparse entries, and residuals are exact rationals. Only +the final evaluation/error comparison uses floating point. + +## 4. Strong-baseline red team + +For 100 queries at `K=20`: + +```text +compile once: + 55 sparse entries + 20 triangular divisions +online compiled proxy: + 9 nonzero series terms per query +strong recurrence: + 9 executed correction steps per query in binary64 +``` + +The online structural counts are comparable. Compilation overhead means the +AMP path does **not** earn a universal single-query floating-point speedup. +At `K=10`, four nonzero terms give about `1.5e-8` absolute error and can be an +economical batch approximation, but the tradeoff depends on tolerance, +initial chart, numeric backend, and number of queries. + +The main earned advantages are instead: + +- avoiding exponential symbolic support; +- obtaining the long-horizon observer without choosing one horizon `N`; +- exact coefficient and residual certificates; +- a reusable compile-once coordinate for many states or parameter sweeps; +- exposing the support geometry and failure boundary. + +## 5. Negative chart control + +At `y_0=0`, the asymptotic series is outside its safe region. The errors are + +```text +K=14: about 0.056 +K=20: more than 4.0 +``` + +Increasing observer order makes the answer worse. A finite AMP truncation is +therefore not a globally convergent numerical method. The compiler must +carry a chart/domain certificate or use residual-driven adaptation; order +alone is not safety. + +## 6. Replay + +The executable tests independently verify: + +- every substitution-matrix coefficient against SymPy series composition; +- nilpotence at the finite observer; +- the exact eigenrelation and first omitted residual; +- expanded support counts against explicit small iterates; +- convergence against the strong logarithmic recurrence; +- sparse/dense cost separation; +- failure outside the asymptotic chart; +- typed refusal for invalid or cancellation-prone frozen tasks. diff --git a/sonnet/amp-polynomial-matrix-compiler/03-disposition.md b/sonnet/amp-polynomial-matrix-compiler/03-disposition.md new file mode 100644 index 00000000..6fea4614 --- /dev/null +++ b/sonnet/amp-polynomial-matrix-compiler/03-disposition.md @@ -0,0 +1,87 @@ +# Disposition: selective algorithmic simplification + +Status: issue #152 result. + +## 1. Verdict by layer + +| Layer | Disposition | Reason | +|---|---|---| +| AMP polynomial-like basis | **EXPAND** | compresses degree-`d^N` expanded support into a fixed observer ray and exposes the correct asymptotic coordinate | +| AMP matrix-like transport | **EXPAND** | composition becomes an exact sparse nilpotent matrix and the conjugacy becomes a triangular linear solve | +| Exact certificate/replay | **EXPAND** | rational coefficients, finite eigenrelation, and first omitted residual replay independently | +| Single-query floating-point speed | **NARROW** | strong logarithmic recurrence stops early and can equal or beat compiled evaluation | +| Global numerical method | **STOP outside chart** | higher truncation can diverge near the non-asymptotic region | +| Generic interacting dynamics | **OPEN** | one power-dominant scalar family is not a general AMP solver | + +The overall disposition is + +\[ +\boxed{\texttt{EXPAND-NARROW}}. +\] + +Both proposed layers earned real algorithmic roles, but only for specified +observers and charts. + +## 2. What was actually simplified + +The representation changes the algorithmic structure: + +\[ +\text{nonlinear repeated state map} +\longrightarrow +\text{sparse linear composition operator} +\longrightarrow +\text{finite eigen-coordinate} +\longrightarrow +\text{one long-horizon observer evaluation}. +\] + +The polynomial-like side supplies the task-adapted dictionary. The +matrix-like side supplies the reusable action and solve. Neither layer earns +the result alone: + +- without the exponential ray, a generic ordinary polynomial/Taylor basis + misses the sparse scale transport; +- without the matrix, the basis remains only a vocabulary and does not compile + iteration into an eigenproblem. + +This is the first exact example in the AMP line where the two sides form one +algorithm rather than two analogies. + +## 3. What remains classical and what is programme-specific + +The Böttcher conjugacy, escape-rate function, and Koopman composition operator +are classical. No priority or replacement claim is earned. + +The Process Geometry contribution under test is the selection principle: + +> start from arithmetic process ranks, choose the chart in which the dominant +> rank is affine, complete the lower-rank residual along its generated support, +> and compile the resulting observer transport. + +In this example that principle independently selects the classical useful +coordinate and produces an exact sparse implementation. More examples are +needed before calling the selection principle general. + +## 4. Next gate + +Do not enlarge the runtime generically. The next discriminating tasks are: + +1. add a mixed logarithmic degree `q^k y^n` so that genuine P/A bracket terms, + not only the `n=0` ray, are required; +2. test a two-variable coupled map where the anchor kernel and cross-variable + support become task-visible; +3. compare automatic AMP dictionary generation against a generic monomial + Carleman dictionary under the same residual and storage budget; +4. require a chart-switch certificate when the asymptotic expansion fails. + +Only transfer to a coupled or mixed-log workload would justify extracting a +reusable engineering abstraction. + +```text +Mathematical Core: unchanged +Research Programme: one T1 positive AMP algorithm calibration +Engineering Architecture: Sonnet-local sparse compiler; no extraction yet +Theory Map: no promoted node +Experimental/Public API: none +``` diff --git a/sonnet/amp-polynomial-matrix-compiler/README.md b/sonnet/amp-polynomial-matrix-compiler/README.md new file mode 100644 index 00000000..fe16a738 --- /dev/null +++ b/sonnet/amp-polynomial-matrix-compiler/README.md @@ -0,0 +1,30 @@ +# AMP polynomial/matrix compiler + +This research-local Sonnet answers issue +[#152](https://github.com/mountain/process-geometry/issues/152). + +Read in order: + +1. [`00-problem-frontier.md`](00-problem-frontier.md) freezes the observer, + baselines, metrics, and claim ceiling. +2. [`01-sparse-compiler-theorems.md`](01-sparse-compiler-theorems.md) derives + the polynomial-like coordinate and matrix-like transport. +3. [`02-benchmark-results.md`](02-benchmark-results.md) records exact support, + sparse cost, numerical error, and the strong-baseline red team. +4. [`03-disposition.md`](03-disposition.md) states where algorithmic + simplification was and was not earned. + +The executable certificate is +[`amp_escape_compiler.py`](amp_escape_compiler.py); its independent tests are +[`test_amp_polynomial_matrix_compiler.py`](../../tests/research/test_amp_polynomial_matrix_compiler.py). + +Current result: + +```text +polynomial-like basis: EXPAND for symbolic support and observer coordinates +matrix-like transport: EXPAND for exact sparse compilation and replay +generic numerical acceleration: NARROW / task and tolerance dependent +overall issue disposition: EXPAND-NARROW +``` + +No Mathematical Core, architecture, dependency, or Public API change follows. diff --git a/sonnet/amp-polynomial-matrix-compiler/amp_escape_compiler.py b/sonnet/amp-polynomial-matrix-compiler/amp_escape_compiler.py new file mode 100644 index 00000000..25864552 --- /dev/null +++ b/sonnet/amp-polynomial-matrix-compiler/amp_escape_compiler.py @@ -0,0 +1,417 @@ +"""AMP polynomial/matrix compiler for a power-dominant iteration. + +For ``f(x)=x**d+t`` on a positive large-``x`` chart, put ``y=log(x)`` and +``q=exp(-y)``. The induced map is + + F(y) = d*y + log(1+t*q**d), + g(q) = q**d / (1+t*q**d). + +On the finite observer basis ``q, ..., q**K``, composition by ``g`` is an +exact sparse matrix ``C``. The truncated escape/Bottcher coordinate + + H_K(y) = y + sum(h[k] q**k) + +is obtained from the *linear* equation ``(d I - C) h = u``. This is the +research-local matrix-like realization of an AMP polynomial-like chart. + +The implementation deliberately uses exact ``Fraction`` arithmetic and keeps +symbolic expansion, strong numerical recurrence, and compiled evaluation as +separate cost baselines. It is not a public Koopman or Bottcher API. +""" + +from __future__ import annotations + +from dataclasses import asdict, dataclass +from fractions import Fraction +from math import comb, exp, log1p + + +def _fraction(value: int | Fraction) -> Fraction: + return value if isinstance(value, Fraction) else Fraction(value) + + +@dataclass(frozen=True) +class SparseSubstitutionMatrix: + """Truncated matrix of ``q**k -> g(q)**k`` with one-based degrees.""" + + order: int + degree: int + interaction: Fraction + entries: tuple[tuple[int, int, Fraction], ...] + + @property + def nonzero_count(self) -> int: + return len(self.entries) + + @property + def dense_entry_count(self) -> int: + return self.order * self.order + + @property + def nilpotence_index_bound(self) -> int: + """Smallest ``r`` guaranteed to satisfy ``C**r=0``.""" + + index = 0 + leading_degree = 1 + while leading_degree <= self.order: + leading_degree *= self.degree + index += 1 + return index + + def apply(self, coefficients: tuple[Fraction, ...]) -> tuple[Fraction, ...]: + if len(coefficients) != self.order: + raise ValueError("coefficient vector has the wrong observer order") + result = [Fraction(0) for _ in range(self.order)] + for row, column, value in self.entries: + result[row - 1] += value * coefficients[column - 1] + return tuple(result) + + +def substitution_coefficient( + row: int, + column: int, + degree: int, + interaction: int | Fraction, +) -> Fraction: + r"""Return ``[q**row] (q**d/(1+t*q**d))**column`` exactly.""" + + if row < 1 or column < 1: + raise ValueError("matrix degrees are one-based positive integers") + if degree < 2: + raise ValueError("the power-dominant gate requires degree >= 2") + + if row % degree or row < degree * column: + return Fraction(0) + tail_index = row // degree - column + t = _fraction(interaction) + return ( + Fraction((-1) ** tail_index * comb(column + tail_index - 1, tail_index)) + * t**tail_index + ) + + +def build_substitution_matrix( + degree: int, + interaction: int | Fraction, + order: int, +) -> SparseSubstitutionMatrix: + if degree < 2: + raise ValueError("degree must be at least two") + if order < 1: + raise ValueError("observer order must be positive") + t = _fraction(interaction) + if t <= 0: + raise ValueError("the frozen benchmark uses a positive interaction") + + entries: list[tuple[int, int, Fraction]] = [] + for column in range(1, order + 1): + tail_index = 0 + while degree * (column + tail_index) <= order: + row = degree * (column + tail_index) + value = substitution_coefficient( + row, + column, + degree, + t, + ) + if value: + entries.append((row, column, value)) + tail_index += 1 + + return SparseSubstitutionMatrix( + order=order, + degree=degree, + interaction=t, + entries=tuple(entries), + ) + + +def interaction_log_coefficients( + degree: int, + interaction: int | Fraction, + order: int, +) -> tuple[Fraction, ...]: + r"""Coefficients of ``log(1+t*q**degree)`` through ``q**order``.""" + + if degree < 2 or order < 1: + raise ValueError("invalid degree or observer order") + t = _fraction(interaction) + result = [Fraction(0) for _ in range(order)] + power = 1 + while degree * power <= order: + result[degree * power - 1] = Fraction((-1) ** (power + 1), power) * t**power + power += 1 + return tuple(result) + + +@dataclass(frozen=True) +class EscapeCoordinate: + """Finite AMP polynomial-like coordinate ``y + sum h_k exp(-k y)``.""" + + degree: int + interaction: Fraction + order: int + coefficients: tuple[Fraction, ...] + + @property + def nonzero_term_count(self) -> int: + return sum(value != 0 for value in self.coefficients) + + def evaluate(self, log_state: float) -> float: + q = exp(-log_state) + power = q + correction = 0.0 + for coefficient in self.coefficients: + correction += float(coefficient) * power + power *= q + return log_state + correction + + +@dataclass(frozen=True) +class ResidualTerm: + degree: int + coefficient: Fraction + + +@dataclass(frozen=True) +class CompilationCost: + observer_order: int + dense_matrix_entries: int + sparse_matrix_entries: int + coordinate_terms: int + triangular_divisions: int + + +@dataclass(frozen=True) +class EscapeCompilationCertificate: + matrix: SparseSubstitutionMatrix + source: tuple[Fraction, ...] + coordinate: EscapeCoordinate + first_omitted_residual: ResidualTerm | None + cost: CompilationCost + + def replay_eigenrelation(self) -> bool: + transported = self.matrix.apply(self.coordinate.coefficients) + d = Fraction(self.coordinate.degree) + return all( + self.source[index] + transported[index] + == d * self.coordinate.coefficients[index] + for index in range(self.coordinate.order) + ) + + +def _residual_terms( + coordinate: EscapeCoordinate, + maximum_degree: int, +) -> tuple[ResidualTerm, ...]: + d = coordinate.degree + t = coordinate.interaction + order = coordinate.order + result: list[ResidualTerm] = [] + + for row in range(1, maximum_degree + 1): + source = Fraction(0) + if row % d == 0: + power = row // d + source = Fraction((-1) ** (power + 1), power) * t**power + + transported = sum( + substitution_coefficient(row, column, d, t) + * coordinate.coefficients[column - 1] + for column in range(1, order + 1) + ) + target = ( + d * coordinate.coefficients[row - 1] + if row <= order + else Fraction(0) + ) + coefficient = source + transported - target + if coefficient: + result.append(ResidualTerm(row, coefficient)) + + return tuple(result) + + +def compile_escape_coordinate( + degree: int, + interaction: int | Fraction, + order: int, +) -> EscapeCompilationCertificate: + """Solve the finite AMP Koopman eigenproblem by triangular transport.""" + + matrix = build_substitution_matrix(degree, interaction, order) + source = interaction_log_coefficients(degree, interaction, order) + + entries_by_row: dict[int, list[tuple[int, Fraction]]] = {} + for row, column, value in matrix.entries: + entries_by_row.setdefault(row, []).append((column, value)) + + coefficients = [Fraction(0) for _ in range(order)] + for row in range(1, order + 1): + transported = sum( + value * coefficients[column - 1] + for column, value in entries_by_row.get(row, ()) + ) + coefficients[row - 1] = (source[row - 1] + transported) / degree + + coordinate = EscapeCoordinate( + degree=degree, + interaction=_fraction(interaction), + order=order, + coefficients=tuple(coefficients), + ) + residuals = _residual_terms(coordinate, maximum_degree=2 * order + degree) + omitted = tuple(term for term in residuals if term.degree > order) + + certificate = EscapeCompilationCertificate( + matrix=matrix, + source=source, + coordinate=coordinate, + first_omitted_residual=omitted[0] if omitted else None, + cost=CompilationCost( + observer_order=order, + dense_matrix_entries=matrix.dense_entry_count, + sparse_matrix_entries=matrix.nonzero_count, + coordinate_terms=coordinate.nonzero_term_count, + triangular_divisions=order, + ), + ) + if not certificate.replay_eigenrelation(): + raise AssertionError("compiled coordinate failed exact replay") + return certificate + + +@dataclass(frozen=True) +class DirectIterationResult: + value: float + executed_steps: int + + +def direct_normalized_log_iteration_result( + degree: int, + interaction: int | Fraction, + initial_log_state: float, + iterations: int, +) -> DirectIterationResult: + """Strong O(N) baseline for ``d**(-N) log(f**N(x))``. + + The normalized value is accumulated without constructing ``d**N`` and the + loop stops safely once the interaction correction underflows to zero. + """ + + if degree < 2 or iterations < 0: + raise ValueError("invalid degree or iteration count") + t = float(_fraction(interaction)) + if t <= 0: + raise ValueError("the frozen benchmark uses a positive interaction") + + current = float(initial_log_state) + normalized = current + inverse_scale = 1.0 + executed_steps = 0 + for _ in range(iterations): + interaction_argument = t * exp(-degree * current) + correction = log1p(interaction_argument) + inverse_scale /= degree + normalized += correction * inverse_scale + current = degree * current + correction + executed_steps += 1 + if correction == 0.0: + break + return DirectIterationResult(normalized, executed_steps) + + +def direct_normalized_log_iteration( + degree: int, + interaction: int | Fraction, + initial_log_state: float, + iterations: int, +) -> float: + return direct_normalized_log_iteration_result( + degree, + interaction, + initial_log_state, + iterations, + ).value + + +def expanded_symbolic_term_count( + degree: int, + interaction: int | Fraction, + iterations: int, +) -> int: + """Exact support count for fully expanded positive ``(x**d+t)`` iterates.""" + + if degree < 2 or iterations < 0: + raise ValueError("invalid degree or iteration count") + t = _fraction(interaction) + if iterations == 0 or t == 0: + return 1 + if t < 0: + raise ValueError("negative interactions may cancel expanded terms") + return degree ** (iterations - 1) + 1 + + +@dataclass(frozen=True) +class BenchmarkReport: + degree: int + interaction: str + observer_order: int + horizon: int + queries: int + initial_log_state: float + direct_normalized_value: float + compiled_normalized_value: float + absolute_error: float + expanded_symbolic_terms: int + direct_recurrence_steps: int + compiled_online_series_terms: int + compilation_cost: dict[str, int] + first_omitted_residual: tuple[int, str] | None + + +def benchmark_report( + *, + degree: int = 2, + interaction: int | Fraction = 1, + observer_order: int = 20, + horizon: int = 100, + queries: int = 1, + initial_log_state: float = 1.5, +) -> BenchmarkReport: + if queries < 1: + raise ValueError("query count must be positive") + certificate = compile_escape_coordinate(degree, interaction, observer_order) + direct_result = direct_normalized_log_iteration_result( + degree, + interaction, + initial_log_state, + horizon, + ) + direct = direct_result.value + compiled = certificate.coordinate.evaluate(initial_log_state) + residual = certificate.first_omitted_residual + return BenchmarkReport( + degree=degree, + interaction=str(_fraction(interaction)), + observer_order=observer_order, + horizon=horizon, + queries=queries, + initial_log_state=initial_log_state, + direct_normalized_value=direct, + compiled_normalized_value=compiled, + absolute_error=abs(direct - compiled), + expanded_symbolic_terms=expanded_symbolic_term_count( + degree, + interaction, + horizon, + ), + direct_recurrence_steps=queries * direct_result.executed_steps, + compiled_online_series_terms=( + queries * certificate.coordinate.nonzero_term_count + ), + compilation_cost=asdict(certificate.cost), + first_omitted_residual=( + (residual.degree, str(residual.coefficient)) if residual else None + ), + ) diff --git a/tests/research/test_amp_polynomial_matrix_compiler.py b/tests/research/test_amp_polynomial_matrix_compiler.py new file mode 100644 index 00000000..80678aad --- /dev/null +++ b/tests/research/test_amp_polynomial_matrix_compiler.py @@ -0,0 +1,180 @@ +"""Exact and numerical audit of the AMP polynomial/matrix compiler.""" + +from __future__ import annotations + +from fractions import Fraction +import importlib.util +from pathlib import Path +import sys + +import sympy as sp + + +MODULE_PATH = ( + Path(__file__).parents[2] + / "sonnet/amp-polynomial-matrix-compiler/amp_escape_compiler.py" +) +SPEC = importlib.util.spec_from_file_location("amp_escape_compiler", MODULE_PATH) +assert SPEC and SPEC.loader +module = importlib.util.module_from_spec(SPEC) +sys.modules[SPEC.name] = module +SPEC.loader.exec_module(module) + + +def test_sparse_matrix_is_exact_substitution_on_the_amp_ray_basis(): + q = sp.symbols("q") + degree = 3 + interaction = sp.Rational(2, 3) + order = 18 + g = q**degree / (1 + interaction * q**degree) + + matrix = module.build_substitution_matrix( + degree, + Fraction(2, 3), + order, + ) + entries = {(row, column): value for row, column, value in matrix.entries} + + for column in range(1, order + 1): + expanded = sp.series(g**column, q, 0, order + 1).removeO().expand() + for row in range(1, order + 1): + expected = expanded.coeff(q, row) + actual = entries.get((row, column), Fraction(0)) + assert expected == sp.Rational(actual.numerator, actual.denominator) + + +def test_substitution_matrix_is_nilpotent_at_a_finite_observer(): + matrix = module.build_substitution_matrix(2, 1, 20) + vector = (Fraction(1),) + (Fraction(0),) * 19 + + for _ in range(matrix.nilpotence_index_bound - 1): + vector = matrix.apply(vector) + assert any(vector) + + vector = matrix.apply(vector) + assert not any(vector) + assert matrix.nilpotence_index_bound == 5 + + +def test_compiled_coordinate_solves_the_exact_finite_eigenproblem(): + certificate = module.compile_escape_coordinate(2, 1, 10) + nonzero = { + degree: coefficient + for degree, coefficient in enumerate( + certificate.coordinate.coefficients, + start=1, + ) + if coefficient + } + + assert nonzero == { + 2: Fraction(1, 2), + 6: Fraction(-1, 3), + 8: Fraction(5, 8), + 10: Fraction(-9, 10), + } + assert certificate.replay_eigenrelation() + assert certificate.first_omitted_residual == module.ResidualTerm( + degree=12, + coefficient=Fraction(2), + ) + + +def test_polynomial_like_support_avoids_expanded_iterate_growth(): + x = sp.symbols("x") + + for degree, maximum_iteration in ((2, 5), (3, 3)): + iterate = x + for iteration in range(maximum_iteration + 1): + actual_terms = len(sp.Poly(sp.expand(iterate), x).terms()) + expected_terms = module.expanded_symbolic_term_count( + degree, + 1, + iteration, + ) + assert actual_terms == expected_terms + iterate = sp.expand(iterate**degree + 1) + + assert module.expanded_symbolic_term_count(2, 1, 100) == 2**99 + 1 + + +def test_compiled_coordinate_converges_to_the_strong_log_recurrence_baseline(): + initial_log_state = 1.5 + direct = module.direct_normalized_log_iteration( + 2, + 1, + initial_log_state, + 200, + ) + + errors = [] + for order in (6, 10, 14, 20): + coordinate = module.compile_escape_coordinate(2, 1, order).coordinate + errors.append(abs(coordinate.evaluate(initial_log_state) - direct)) + + assert errors[0] < 4e-6 + assert errors[1] < 2e-8 + assert errors[2] < 5e-11 + assert errors[3] < 5e-13 + assert errors == sorted(errors, reverse=True) + + +def test_sparse_cost_is_reported_without_hiding_the_strong_baseline(): + report = module.benchmark_report( + degree=2, + interaction=1, + observer_order=20, + horizon=100, + queries=100, + initial_log_state=1.5, + ) + + assert report.compilation_cost == { + "observer_order": 20, + "dense_matrix_entries": 400, + "sparse_matrix_entries": 55, + "coordinate_terms": 9, + "triangular_divisions": 20, + } + assert report.expanded_symbolic_terms == 2**99 + 1 + assert report.compiled_online_series_terms == 900 + # The strong numerical baseline detects underflow of the correction and + # stops early; the compiler therefore does not receive a false O(100) + # per-query advantage on this floating-point task. + assert report.direct_recurrence_steps < 2_000 + assert report.absolute_error < 5e-13 + assert report.first_omitted_residual == (22, "144/11") + + +def test_asymptotic_chart_failure_is_not_hidden_by_higher_order(): + initial_log_state = 0.0 + direct = module.direct_normalized_log_iteration(2, 1, initial_log_state, 200) + order_14 = module.compile_escape_coordinate(2, 1, 14).coordinate.evaluate( + initial_log_state + ) + order_20 = module.compile_escape_coordinate(2, 1, 20).coordinate.evaluate( + initial_log_state + ) + + error_14 = abs(order_14 - direct) + error_20 = abs(order_20 - direct) + assert error_14 < 0.06 + assert error_20 > 4.0 + assert error_20 > error_14 + + +def test_invalid_or_cancellation_prone_tasks_fail_closed(): + for args in ((1, 1, 10), (2, 0, 10), (2, -1, 10), (2, 1, 0)): + try: + module.build_substitution_matrix(*args) + except ValueError: + pass + else: # pragma: no cover + raise AssertionError(f"invalid compiler task was accepted: {args}") + + try: + module.expanded_symbolic_term_count(2, -1, 5) + except ValueError: + pass + else: # pragma: no cover + raise AssertionError("a cancellation-prone term count was accepted")