mirror of
https://github.com/allaunthefox/Research-Stack.git
synced 2026-07-31 03:05:21 +00:00
allaun
475f6319ea
This squashes all local history (768 commits) onto the scrubbed PR #90 baseline. Individual commits were lost during filter-repo corruption; the working tree content is preserved intact. Build: N/A (working tree state only)
238 lines
8 KiB
Python
238 lines
8 KiB
Python
#!/usr/bin/env python3
|
||
"""
|
||
QR → Braid Bridge: Householder QR decomposition mapped to 8×8 braid crossings.
|
||
|
||
Each Householder reflector H = I - τ·v·vᵀ is structurally isomorphic to a
|
||
braid crossing operator on the 8-strand braid. The product Q = H₁·…·Hₖ
|
||
is the braid word; R is the upper-triangular "residual" after all crossings
|
||
are applied.
|
||
|
||
For an 8×8 matrix (our Matrix8 type), QR decomposition applies exactly
|
||
8 Householder reflectors — one per strand crossing in the braid.
|
||
|
||
Usage:
|
||
python3 qr_braid_bridge.py
|
||
"""
|
||
|
||
import json
|
||
import hashlib
|
||
import math
|
||
import random
|
||
from datetime import datetime, timezone
|
||
from dataclasses import dataclass
|
||
|
||
EPSILON = 1e-14
|
||
|
||
|
||
@dataclass
|
||
class HouseholderReflector:
|
||
"""A single Householder reflector = one braid crossing."""
|
||
k: int # column being zeroed (strand index)
|
||
v: list[float] # reflector vector v (normalized)
|
||
tau: float # scaling factor τ
|
||
|
||
|
||
def householder_vector(x: list[float]) -> tuple[list[float], float]:
|
||
"""Compute Householder reflector vector v and τ for vector x.
|
||
H = I - τ·v·vᵀ such that H·x = μ·e₀."""
|
||
n = len(x)
|
||
alpha = x[0]
|
||
norm_x = math.sqrt(sum(xi * xi for xi in x))
|
||
mu = -math.copysign(norm_x, alpha)
|
||
|
||
if abs(mu - alpha) < EPSILON:
|
||
return [0.0] * n, 0.0 # no reflection needed
|
||
|
||
v = [0.0] * n
|
||
v[0] = 1.0
|
||
for i in range(1, n):
|
||
v[i] = x[i] / (alpha - mu)
|
||
tau = (mu - alpha) / mu
|
||
return v, tau
|
||
|
||
|
||
def apply_householder(A: list[list[float]], k: int, v: list[float], tau: float):
|
||
"""Apply Householder reflector H = I - τ·v·vᵀ to trailing submatrix A[k:, k:]."""
|
||
m = len(A)
|
||
n = len(A[0])
|
||
|
||
# Apply to rows k..m-1, columns k..n-1
|
||
# H·A = A - τ·v·(vᵀ·A)
|
||
for j in range(k, n):
|
||
dot = sum(v[i] * A[k + i][j] for i in range(m - k))
|
||
for i in range(m - k):
|
||
A[k + i][j] -= tau * v[i] * dot
|
||
|
||
|
||
def qr_decomposition(A: list[list[float]]) -> tuple[list[list[float]], list[list[float]], list[HouseholderReflector]]:
|
||
"""QR decomposition via Householder reflections.
|
||
Returns Q, R, and the list of Householder reflectors (braid crossings)."""
|
||
m = len(A)
|
||
n = len(A[0])
|
||
k = min(m, n)
|
||
|
||
R = [row[:] for row in A] # copy
|
||
reflectors = []
|
||
|
||
for col in range(k):
|
||
# Extract column x = R[col:, col]
|
||
x = [R[i][col] for i in range(col, m)]
|
||
v_raw, tau = householder_vector(x)
|
||
|
||
if abs(tau) < EPSILON:
|
||
continue
|
||
|
||
# v has length m-col; pad with zeros at the top for the full matrix
|
||
v_full = [0.0] * col + v_raw
|
||
|
||
reflectors.append(HouseholderReflector(k=col, v=v_full, tau=tau))
|
||
|
||
# Apply reflector to R
|
||
for j in range(col, n):
|
||
dot = sum(v_full[i] * R[i][j] for i in range(col, m))
|
||
for i in range(col, m):
|
||
R[i][j] -= tau * v_full[i] * dot
|
||
|
||
# Construct Q = H₁·H₂·...·Hₖ
|
||
Q = [[float(i == j) for j in range(m)] for i in range(m)]
|
||
for ref in reversed(reflectors):
|
||
for j in range(m):
|
||
dot = sum(ref.v[i] * Q[i][j] for i in range(m))
|
||
for i in range(m):
|
||
Q[i][j] -= ref.tau * ref.v[i] * dot
|
||
|
||
# Zero out near-zero elements
|
||
for i in range(m):
|
||
for j in range(n):
|
||
if abs(R[i][j]) < EPSILON:
|
||
R[i][j] = 0.0
|
||
if abs(Q[i][j]) < EPSILON:
|
||
Q[i][j] = 0.0
|
||
|
||
return Q, R, reflectors
|
||
|
||
|
||
def matrix8_example() -> list[list[float]]:
|
||
"""Our canonical test matrix from AdjugateMatrix.lean: identity.lean."""
|
||
return [[1, 0, 0, 0, 0, 0, 0, 0],
|
||
[0, 1, 0, 0, 0, 0, 0, 0],
|
||
[0, 0, 1, 0, 0, 0, 0, 0],
|
||
[0, 0, 0, 1, 0, 0, 0, 0],
|
||
[0, 0, 0, 0, 1, 0, 0, 0],
|
||
[0, 0, 0, 0, 0, 1, 0, 0],
|
||
[0, 0, 0, 0, 0, 0, 1, 0],
|
||
[0, 0, 0, 0, 0, 0, 0, 1]]
|
||
|
||
|
||
def matrix8_diag3() -> list[list[float]]:
|
||
"""The Q16_16 counterexample from AdjugateMatrix.lean: diag(3,1,...,1)."""
|
||
m = [[0] * 8 for _ in range(8)]
|
||
for i in range(8):
|
||
m[i][i] = 3.0 if i == 0 else 1.0
|
||
return m
|
||
|
||
|
||
def matrix8_random() -> list[list[float]]:
|
||
"""Random 8×8 matrix for testing."""
|
||
return [[random.uniform(-1, 1) for _ in range(8)] for _ in range(8)]
|
||
|
||
|
||
def mat_mul(A: list[list[float]], B: list[list[float]]) -> list[list[float]]:
|
||
"""Standard matrix multiply."""
|
||
m, n, p = len(A), len(A[0]), len(B[0])
|
||
return [[sum(A[i][k] * B[k][j] for k in range(n)) for j in range(p)] for i in range(m)]
|
||
|
||
|
||
def frobenius_error(A: list[list[float]], B: list[list[float]]) -> float:
|
||
"""Frobenius norm of A - B."""
|
||
m, n = len(A), len(A[0])
|
||
return math.sqrt(sum((A[i][j] - B[i][j])**2 for i in range(m) for j in range(n)))
|
||
|
||
|
||
def sidon_weight(ref: HouseholderReflector) -> int:
|
||
"""Map a Householder reflector (column k) to a Sidon weight.
|
||
For 8×8, columns 0..7 map to Sidon addresses {1,2,4,8,16,32,64,128}."""
|
||
return 1 << ref.k
|
||
|
||
|
||
def main():
|
||
print("=" * 60)
|
||
print("QR → Braid Bridge: Householder = Crossing, WY = Strand Block")
|
||
print("=" * 60)
|
||
|
||
test_matrices = [
|
||
("Identity (I₈)", matrix8_example()),
|
||
("diag(3,1,1,1,1,1,1,1)", matrix8_diag3()),
|
||
("Random 8×8 #1", matrix8_random()),
|
||
("Random 8×8 #2", matrix8_random()),
|
||
]
|
||
|
||
results = []
|
||
for name, A in test_matrices:
|
||
print(f"\n{'─'*60}")
|
||
print(f"Matrix: {name}")
|
||
|
||
Q, R, reflectors = qr_decomposition(A)
|
||
A_reconstructed = mat_mul(Q, R)
|
||
error = frobenius_error(A, A_reconstructed)
|
||
n_reflectors = len(reflectors)
|
||
|
||
# Compute Sidon weights
|
||
sidon_weights = [sidon_weight(ref) for ref in reflectors]
|
||
total_sidon = sum(sidon_weights)
|
||
sidon_max = 128 # for 8×8, max Sidon label is 2^7 = 128
|
||
|
||
print(f" QR error (||QR - A||_F): {error:.2e}")
|
||
print(f" Householder reflectors: {n_reflectors} (braid crossings)")
|
||
print(f" Sidon weights: {sidon_weights}")
|
||
print(f" Total Sidon weight: {total_sidon} / {sidon_max} (sidon_slack={sidon_max - total_sidon})")
|
||
|
||
for i, ref in enumerate(reflectors):
|
||
nonzero = sum(1 for vi in ref.v if abs(vi) > EPSILON)
|
||
print(f" H[{i}] (strand {ref.k}): tau={ref.tau:.4f}, {nonzero} nonzero components")
|
||
|
||
results.append({
|
||
"name": name,
|
||
"qr_error": round(error, 10),
|
||
"n_reflectors": n_reflectors,
|
||
"sidon_weights": sidon_weights,
|
||
"total_sidon_weight": total_sidon,
|
||
"sidon_slack": sidon_max - total_sidon,
|
||
})
|
||
|
||
# Build receipt
|
||
receipt = {
|
||
"schema": "rrc_qr_braid_bridge_v1",
|
||
"claim_boundary": "householder_qr_on_8x8;each_reflector_is_braid_crossing",
|
||
"description": (
|
||
"QR decomposition of an 8×8 matrix requires at most 8 Householder "
|
||
"reflectors, one per column. Each reflector maps to a braid crossing "
|
||
"with Sidon weight 2^k for column k. The total Sidon weight is the "
|
||
"sum of active crossing addresses; sidon_slack measures unused "
|
||
"crossing capacity."
|
||
),
|
||
"tests": results,
|
||
"summary": {
|
||
"total_tests": len(results),
|
||
"max_error": max(r["qr_error"] for r in results),
|
||
"max_reflectors": max(r["n_reflectors"] for r in results),
|
||
"braid_word_length": f"≤8 crossings per 8×8 decomposition",
|
||
"sidon_budget": 128,
|
||
},
|
||
"computed_at": datetime.now(timezone.utc).isoformat(),
|
||
}
|
||
|
||
canonical = json.dumps(receipt, sort_keys=True, separators=(",", ":"))
|
||
receipt["receipt_sha256"] = hashlib.sha256(canonical.encode()).hexdigest()
|
||
|
||
path = "qr_braid_bridge_receipt.json"
|
||
with open(path, "w") as f:
|
||
json.dump(receipt, f, indent=2, sort_keys=True)
|
||
|
||
print(f"\n{'='*60}")
|
||
print(f"Receipt: {path}")
|
||
print(f"SHA256: {receipt['receipt_sha256']}")
|
||
|
||
|
||
if __name__ == "__main__":
|
||
main()
|