Create prototype.py

prototype.py

@@ -0,0 +1,348 @@
+"""
+FLEighConduit — Evidence-Emitting Eigendecomposition
+======================================================
+Extends FLEigh with judicial conduit telemetry.
+Read-only observation of the solver's internal states.
+
+Council-ratified specification (11 rounds, 3 AI participants):
+  - Theorem 1: Lens Preservation (shared arithmetic path)
+  - Theorem 2: Dynamic Non-Reconstructibility (friction, settle, order)
+  - Theorem 4: Continuity (static continuous, dynamic piecewise)
+  - Theorem 5: Gauge-Safe Directional Observation (sign canonicalization)
+
+Classes:
+  ConduitPacket  — fixed-shape tensor bundle, batch-first, dimension-agnostic
+  FLEighConduit  — extends FLEigh, emits ConduitPacket
+
+Usage:
+    from geolip_core.linalg.conduit import FLEighConduit
+
+    solver = FLEighConduit()
+    packet = solver(A)          # A: (B, n, n) symmetric
+    packet.eigenvalues          # (B, n)
+    packet.friction             # (B, n) — per-root solver struggle
+    packet.settle               # (B, n) — iterations to convergence
+    packet.extraction_order     # (B, n) — root extraction sequence
+
+    # Standard eigenpairs (identical to FLEigh)
+    evals, evecs = packet.eigenpairs()
+
+License: MIT
+
+Author: AbstractPhil + Claude 4.6 Opus Extended
+Assistants: Gemini Pro, GPT 5.4 Extended Thinking
+"""
+
+import math
+import torch
+import torch.nn as nn
+from torch import Tensor
+from typing import Tuple, Optional
+from dataclasses import dataclass
+
+
+@dataclass
+class ConduitPacket:
+    """Fixed-shape telemetry from FLEighConduit.
+
+    All tensors batch-first, dimension-agnostic.
+    Production packet: scalar-dominant, bounded overhead.
+    Research fields (Mstore, trajectories) populated only when requested.
+    """
+
+    # ── Spectral evidence (static, deterministic) ──
+    eigenvalues: Tensor          # (B, n) sorted ascending
+    eigenvectors: Tensor         # (B, n, n) sign-canonicalized
+    char_coeffs: Tensor          # (B, n) elementary symmetric polys, monic 1 omitted
+
+    # ── Adjudication evidence (dynamic, non-reconstructible) ──
+    friction: Tensor             # (B, n) per-root Σ 1/(|p'(z_t)| + δ)
+    settle: Tensor               # (B, n) iterations to convergence per root
+    extraction_order: Tensor     # (B, n) which root found first (0-indexed)
+    refinement_residual: Tensor  # (B,) ||V^T V - I||_F after Newton-Schulz
+
+    # ── Release fidelity ──
+    # Note: release_residual is owned by SVDConduit layer in full architecture.
+    # Included here for v1 convenience when enc_out matrix M is provided.
+    release_residual: Optional[Tensor] = None  # (B,) ||M - U diag(S) Vt||²
+
+    # ── Research mode (populated only with research=True) ──
+    mstore: Optional[Tensor] = None     # (n+1, B, n, n) FL matrix states
+    z_trajectory: Optional[Tensor] = None  # (B, n, laguerre_iters) root guesses
+    dp_trajectory: Optional[Tensor] = None # (B, n, laguerre_iters) p' at each step
+
+    def eigenpairs(self) -> Tuple[Tensor, Tensor]:
+        """Standard output matching FLEigh contract."""
+        return self.eigenvalues, self.eigenvectors
+
+
+def canonicalize_eigenvectors(V: Tensor) -> Tensor:
+    """Force deterministic sign convention on eigenvector columns.
+
+    For each column (eigenvector), flip sign so the entry with
+    largest absolute value is positive. Resolves the gauge ambiguity
+    that otherwise causes identical matrices to produce different
+    embeddings on S^(n²-1).
+
+    Args:
+        V: (B, n, n) eigenvector matrix (columns are eigenvectors)
+    Returns:
+        V with deterministic signs
+    """
+    # Find index of max absolute value per column
+    max_idx = V.abs().argmax(dim=-2, keepdim=True)  # (B, 1, n)
+    sign = V.gather(-2, max_idx).sign()               # (B, 1, n)
+    return V * sign
+
+
+class FLEighConduit(nn.Module):
+    """Evidence-emitting eigendecomposition.
+
+    Identical arithmetic to FLEigh. Captures telemetry at phase
+    boundaries without altering the numerical path.
+
+    Phases (shared with FLEigh):
+      1. FL characteristic polynomial (fp64, n bmm)
+      2. Laguerre root-finding + Newton polish (with telemetry capture)
+      3. FL adjugate eigenvectors (fp64 Horner + max-col)
+      4. Newton-Schulz orthogonalization (fp32, 2 iters)
+      5. Rayleigh quotient refinement (fp32, 2 bmm)
+
+    Args:
+        laguerre_iters: Root-finding iterations per eigenvalue (default 5)
+        polish_iters: Newton refinement iterations (default 3)
+        ns_iters: Newton-Schulz orthogonalization iterations (default 2)
+        friction_delta: stability constant for friction computation (default 1e-8)
+        settle_threshold: convergence threshold for settle count (default 1e-6)
+        research: if True, populate full trajectory and Mstore fields
+    """
+
+    def __init__(self, laguerre_iters: int = 5, polish_iters: int = 3,
+                 ns_iters: int = 2, friction_delta: float = 1e-8,
+                 settle_threshold: float = 1e-6, research: bool = False):
+        super().__init__()
+        self.laguerre_iters = laguerre_iters
+        self.polish_iters = polish_iters
+        self.ns_iters = ns_iters
+        self.friction_delta = friction_delta
+        self.settle_threshold = settle_threshold
+        self.research = research
+
+    def forward(self, A: Tensor) -> ConduitPacket:
+        """Evidence-emitting eigendecomposition.
+
+        Args:
+            A: (B, n, n) symmetric matrix batch
+
+        Returns:
+            ConduitPacket with eigenpairs + judicial telemetry
+        """
+        B, n, _ = A.shape
+        device = A.device
+
+        # ════════════════════════════════════════════
+        # Phase 1: Faddeev-LeVerrier (fp64)
+        # ════════════════════════════════════════════
+        scale = (torch.linalg.norm(A.reshape(B, -1), dim=-1) / math.sqrt(n)).clamp(min=1e-12)
+        As = A / scale[:, None, None]
+
+        Ad = As.double()
+        eye_d = torch.eye(n, device=device, dtype=torch.float64).unsqueeze(0).expand(B, -1, -1)
+        c = torch.zeros(B, n + 1, device=device, dtype=torch.float64)
+        c[:, n] = 1.0
+        Mstore = torch.zeros(n + 1, B, n, n, device=device, dtype=torch.float64)
+        Mk = torch.zeros(B, n, n, device=device, dtype=torch.float64)
+        for k in range(1, n + 1):
+            Mk = torch.bmm(Ad, Mk) + c[:, n - k + 1, None, None] * eye_d
+            Mstore[k] = Mk
+            c[:, n - k] = -(Ad * Mk).sum((-2, -1)) / k
+
+        # Capture characteristic coefficients (omit monic leading 1)
+        char_coeffs = c[:, :n].float()  # (B, n)
+
+        # ════════════════════════════════════════════
+        # Phase 2: Laguerre + deflation + Newton polish
+        #          WITH TELEMETRY CAPTURE
+        # ════════════════════════════════════════════
+        use_f64 = n > 6
+        dt = torch.float64 if use_f64 else torch.float32
+        cl = c.to(dt).clone().detach()
+        roots = torch.zeros(B, n, device=device, dtype=dt)
+        zi = As.to(dt).diagonal(dim1=-2, dim2=-1).sort(dim=-1).values.detach()
+        zi = zi + torch.linspace(-1e-4, 1e-4, n, device=device, dtype=dt).unsqueeze(0)
+
+        # Telemetry buffers
+        friction = torch.zeros(B, n, device=device, dtype=torch.float32)
+        settle = torch.full((B, n), float(self.laguerre_iters),
+                           device=device, dtype=torch.float32)
+        extraction_order = torch.zeros(B, n, device=device, dtype=torch.float32)
+
+        # Research buffers
+        if self.research:
+            z_traj = torch.zeros(B, n, self.laguerre_iters,
+                                device=device, dtype=torch.float32)
+            dp_traj = torch.zeros(B, n, self.laguerre_iters,
+                                 device=device, dtype=torch.float32)
+
+        for ri in range(n):
+            deg = n - ri
+            z = zi[:, ri]
+
+            for lag_iter in range(self.laguerre_iters):
+                # Horner evaluation
+                pv = cl[:, deg]
+                dp = torch.zeros(B, device=device, dtype=dt)
+                d2 = torch.zeros(B, device=device, dtype=dt)
+                for j in range(deg - 1, -1, -1):
+                    d2 = d2 * z + dp
+                    dp = dp * z + pv
+                    pv = pv * z + cl[:, j]
+
+                # ── Telemetry capture ──
+                dp_abs = dp.abs().float()
+                friction[:, ri] += 1.0 / (dp_abs + self.friction_delta)
+
+                # Settle detection
+                pv_abs = pv.abs().float()
+                just_settled = (pv_abs < self.settle_threshold) & \
+                               (settle[:, ri] == float(self.laguerre_iters))
+                settle[:, ri] = torch.where(just_settled,
+                                           torch.full_like(settle[:, ri], float(lag_iter)),
+                                           settle[:, ri])
+
+                if self.research:
+                    z_traj[:, ri, lag_iter] = z.float()
+                    dp_traj[:, ri, lag_iter] = dp_abs
+
+                # Laguerre step (unchanged arithmetic)
+                ok = pv.abs() > 1e-30
+                ps = torch.where(ok, pv, torch.ones_like(pv))
+                G = torch.where(ok, dp / ps, torch.zeros_like(dp))
+                H = G * G - torch.where(ok, 2.0 * d2 / ps, torch.zeros_like(d2))
+                disc = ((deg - 1.0) * (deg * H - G * G)).clamp(min=0.0)
+                sq = torch.sqrt(disc)
+                gp = G + sq
+                gm = G - sq
+                den = torch.where(gp.abs() >= gm.abs(), gp, gm)
+                dok = den.abs() > 1e-20
+                ds = torch.where(dok, den, torch.ones_like(den))
+                z = z - torch.where(dok, float(deg) / ds, torch.zeros_like(den))
+
+            roots[:, ri] = z
+            extraction_order[:, ri] = float(ri)
+
+            # Synthetic division (deflation)
+            b = cl[:, deg]
+            for j in range(deg - 1, 0, -1):
+                bn = cl[:, j] + z * b
+                cl[:, j] = b
+                b = bn
+            cl[:, 0] = b
+
+        # Newton polish on original polynomial
+        roots = roots.double()
+        for _ in range(self.polish_iters):
+            pv = torch.ones(B, n, device=device, dtype=torch.float64)
+            dp = torch.zeros(B, n, device=device, dtype=torch.float64)
+            for j in range(n - 1, -1, -1):
+                dp = dp * roots + pv
+                pv = pv * roots + c[:, j:j + 1]
+            ok = dp.abs() > 1e-30
+            dps = torch.where(ok, dp, torch.ones_like(dp))
+            roots = roots - torch.where(ok, pv / dps, torch.zeros_like(pv))
+
+        # ════════════════════════════════════════════
+        # Phase 3: FL adjugate eigenvectors (fp64)
+        # ════════════════════════════════════════════
+        lam = roots
+        R = Mstore[1].unsqueeze(1).expand(-1, n, -1, -1).clone()
+        for k in range(2, n + 1):
+            R = R * lam[:, :, None, None] + Mstore[k].unsqueeze(1)
+        cnorms = R.norm(dim=-2)
+        best = cnorms.argmax(dim=-1)
+        idx = best.unsqueeze(-1).unsqueeze(-1).expand(-1, -1, n, 1)
+        vec = R.gather(-1, idx).squeeze(-1)
+        vec = vec / (vec.norm(dim=-1, keepdim=True) + 1e-30)
+        V = vec.float().transpose(-2, -1)
+
+        # ════════════════════════════════════════════
+        # Phase 4: Newton-Schulz orthogonalization
+        # ════════════════════════════════════════════
+        eye_f = torch.eye(n, device=device, dtype=torch.float32).unsqueeze(0).expand(B, -1, -1)
+        Y = torch.bmm(V.transpose(-2, -1), V)
+        X = eye_f.clone()
+        for _ in range(self.ns_iters):
+            T = 3.0 * eye_f - Y
+            X = 0.5 * torch.bmm(X, T)
+            Y = 0.5 * torch.bmm(T, Y)
+        V = torch.bmm(V, X)
+
+        # Telemetry: orthogonality residual after NS
+        VtV = torch.bmm(V.transpose(-2, -1), V)
+        refinement_residual = (VtV - eye_f).pow(2).sum((-2, -1)).sqrt()  # (B,)
+
+        # ════════════════════════════════════════════
+        # Phase 5: Rayleigh quotient refinement
+        # ════════════════════════════════════════════
+        AV = torch.bmm(A, V)
+        evals = (V * AV).sum(dim=-2)
+
+        se, perm = evals.sort(dim=-1)
+        sv = V.gather(-1, perm.unsqueeze(-2).expand_as(V))
+
+        # Reorder telemetry to match sorted eigenvalue order
+        friction_sorted = friction.gather(-1, perm)
+        settle_sorted = settle.gather(-1, perm)
+        # extraction_order stays as-is — it records the ORIGINAL extraction sequence
+
+        # ════════════════════════════════════════════
+        # Gauge canonicalization
+        # ════════════════════════════════════════════
+        sv = canonicalize_eigenvectors(sv)
+
+        # ════════════════════════════════════════════
+        # Build packet
+        # ════════════════════════════════════════════
+        packet = ConduitPacket(
+            eigenvalues=se,
+            eigenvectors=sv,
+            char_coeffs=char_coeffs,
+            friction=friction_sorted,
+            settle=settle_sorted,
+            extraction_order=extraction_order,
+            refinement_residual=refinement_residual,
+        )
+
+        if self.research:
+            packet.mstore = Mstore
+            packet.z_trajectory = z_traj
+            packet.dp_trajectory = dp_traj
+
+        return packet
+
+
+# ── Regression parity test ──
+
+def verify_parity(A: Tensor, atol: float = 1e-5) -> bool:
+    """Verify FLEighConduit produces identical eigenpairs to FLEigh.
+
+    Args:
+        A: (B, n, n) symmetric test matrices
+        atol: absolute tolerance
+
+    Returns:
+        True if eigenpairs match within tolerance
+    """
+    from geolip_core.linalg.eigh import FLEigh
+
+    ref_evals, ref_evecs = FLEigh()(A)
+    packet = FLEighConduit()(A)
+    cond_evals, cond_evecs = packet.eigenpairs()
+
+    evals_match = torch.allclose(ref_evals, cond_evals, atol=atol)
+
+    # Eigenvectors may differ by sign — compare via absolute inner products
+    dots = (ref_evecs * cond_evecs).sum(dim=-2).abs()
+    evecs_match = (dots > 1.0 - atol).all()
+
+    return evals_match and evecs_match