1import json, time
  2import numpy as np
  3
  4
  5def rank_project(X, q):
  6    U, s, Vt = np.linalg.svd(X, full_matrices=False)
  7    return (U[:, :q] * s[:q]) @ Vt[:q, :], U[:, :q], Vt[:q, :].T, s
  8
  9
 10def tangent_project(H, U, V):
 11    UtH = U.T @ H
 12    HV = H @ V
 13    UtHV = U.T @ H @ V
 14    return U @ UtH + HV @ V.T - U @ UtHV @ V.T
 15
 16
 17def normal_project(H, U, V):
 18    return H - tangent_project(H, U, V)
 19
 20
 21def adaptive_step(X, B, b, q, c=1.0, rho=1.0, eps=1e-6):
 22    Y, U, V, s = rank_project(X, q)
 23    r = np.linalg.norm(X-Y) / (np.linalg.norm(X) + 1e-12)
 24    mu = c * (r + eps) ** rho
 25    d = X.size
 26    # Explicit Q is appropriate for this deliberately small numerical check.
 27    Q = np.zeros((d, d))
 28    for j in range(d):
 29        E = np.zeros((X.shape[0], X.shape[1])); E.flat[j] = 1.0
 30        Q[:, j] = (normal_project(E, U, V) + mu * E).ravel()
 31    K = np.block([[Q, B.T], [B, np.zeros((B.shape[0], B.shape[0]))]])
 32    rhs = np.concatenate([(mu * X + normal_project(Y, U, V)).ravel(), b])
 33    sol = np.linalg.solve(K, rhs)
 34    Z = sol[:d].reshape(X.shape)
 35    return Z, Y, r, mu
 36
 37
 38def affine_projection(X, B, b):
 39    # Euclidean projection onto B vec(Z)=b.
 40    BBt = B @ B.T
 41    return (X.ravel() + B.T @ np.linalg.solve(BBt, b - B @ X.ravel())).reshape(X.shape)
 42
 43
 44def verify_math(rng):
 45    m, n, q = 7, 6, 2
 46    X = rng.normal(size=(m, n))
 47    Y, U, V, _ = rank_project(X, q)
 48    H = rng.normal(size=(m, n))
 49    PT = tangent_project(H, U, V); PN = normal_project(H, U, V)
 50    projector_err = max(np.linalg.norm(tangent_project(PT, U, V)-PT),
 51                        np.linalg.norm(normal_project(PN, U, V)-PN))
 52    orth_err = abs(np.sum(PT * PN))
 53    decomp_err = np.linalg.norm(PT + PN - H)
 54    # The residual rule should have precisely the claimed power-law log slope.
 55    rs = np.logspace(-6, -1, 20)
 56    mus = np.array([(rs[i] + 1e-12) ** 0.7 for i in range(len(rs))])
 57    slope = np.polyfit(np.log(rs), np.log(mus), 1)[0]
 58    return {"projector_idempotence_error": float(projector_err),
 59            "projector_orthogonality_error": float(orth_err),
 60            "projector_decomposition_error": float(decomp_err),
 61            "mu_log_slope_rho_0.7": float(slope),
 62            "math_pass": bool(projector_err < 1e-10 and orth_err < 1e-10 and decomp_err < 1e-10 and abs(slope-.7)<1e-3)}
 63
 64
 65def run_method(name, W0, target, B, b, q, steps=100, c=1.0, rho=1.0):
 66    W = W0.copy(); hist = []; t0 = time.perf_counter()
 67    for k in range(steps):
 68        # Same gradient step for every method.
 69        X = W - 0.35 * (W - target)
 70        if name == "alternating":
 71            Y, _, _, _ = rank_project(X, q)
 72            W = affine_projection(Y, B, b)
 73            # Keep the standard alternating scheme's next iterate rank feasible.
 74            W, _, _, _ = rank_project(W, q)
 75            mu = None; r = np.linalg.norm(X-Y)/(np.linalg.norm(X)+1e-12)
 76        else:
 77            use_rho = 0.0 if name == "fixed_mu" else rho
 78            Z, Y, r, mu = adaptive_step(X, B, b, q, c=c, rho=use_rho)
 79            W, _, _, _ = rank_project(Z, q)
 80        rank_res = np.linalg.norm(W - rank_project(W, q)[0])/(np.linalg.norm(W)+1e-12)
 81        aff_res = np.linalg.norm(B @ W.ravel() - b)
 82        err = np.linalg.norm(W-target)/np.linalg.norm(target)
 83        hist.append((err, rank_res, aff_res, r, -1 if mu is None else mu))
 84    return {"final_relative_target_error": float(hist[-1][0]),
 85            "final_affine_residual": float(hist[-1][2]),
 86            "final_rank_residual": float(hist[-1][1]),
 87            "initial_affine_residual": float(hist[0][2]),
 88            "min_error_last_20": float(min(x[0] for x in hist[-20:])),
 89            "time_sec": time.perf_counter()-t0,
 90            "history": np.array(hist).tolist()}
 91
 92
 93def main():
 94    rng = np.random.default_rng(284)
 95    math = verify_math(rng)
 96    m, n, q = 8, 7, 2
 97    # A rank-q target defines a consistent affine calibration constraint.
 98    target = rng.normal(size=(m,q)) @ rng.normal(size=(q,n))
 99    d = m*n
100    # Full-row-rank random linear calibration observations.
101    B0 = rng.normal(size=(12,d)); Qb, _ = np.linalg.qr(B0.T); B = Qb[:, :12].T
102    b = B @ target.ravel()
103    W0 = target + 1.8*rng.normal(size=(m,n))
104    configs = {"fixed_mu": (1.0, 0.0), "adaptive_rho05": (1.0, 0.5), "adaptive_rho1": (1.0, 1.0), "adaptive_c10_rho1": (10.0, 1.0), "alternating": (0.0, 0.0)}
105    results = {name: run_method(name, W0, target, B, b, q, c=c, rho=rho) for name, (c, rho) in configs.items()}
106    # Strip histories from summary but save them for reproducibility.
107    out = {"math": math, "benchmark": results,
108           "setup": {"shape": [m,n], "rank": q, "constraints": B.shape[0], "steps": 100,
109                     "learning_rate": .35, "seed": 284}}
110    with open("results.json", "w") as f: json.dump(out, f, indent=2)
111    print(json.dumps({"math": math, "benchmark_summary":
112        {k:{x:v for x,v in val.items() if x != "history"} for k,val in results.items()}}, indent=2))
113
114if __name__ == "__main__": main()