fig_theorem31.py

"""Figures for the Theorem 3.1 audit (Claim 2)."""

import csv
import json

import numpy as np
import plotly.graph_objects as go
from plotly.subplots import make_subplots

RNG = np.random.default_rng(0)
LOG2PI = np.log(2 * np.pi)


def tcr(F, eps=0.5):
    k, d = F.shape
    Fn = F / np.linalg.norm(F, axis=1, keepdims=True)
    return float(-0.5 * np.linalg.slogdet(
        np.eye(k) + (d / (k * eps ** 2)) * (Fn @ Fn.T))[1])


def gauss_tc(cov):
    return float(0.5 * np.log(np.diag(cov)).sum() - 0.5 * np.linalg.slogdet(cov)[1])


k, d = 8, 128
rhos = np.linspace(0.0, 0.95, 40)
tcs, divs = [], []
for rho in rhos:
    cov = (1 - rho) * np.eye(k) + rho * np.ones((k, k))
    F = np.linalg.cholesky(cov) @ RNG.normal(size=(k, d))
    tcs.append(gauss_tc(np.cov(F)))
    divs.append(tcr(F))
tcs, divs = np.array(tcs), np.array(divs)

# panel B: MSE vs Gaussian ELBO (exact) and vs Laplace log-likelihood (control)
n, dd, s2, b = 20000, 6, 0.7, 0.6
x = RNG.normal(size=(n, dd))
scales = np.linspace(0.05, 2.0, 25)
mses, gll, lll = [], [], []
for s in scales:
    xh = x + s * RNG.normal(size=(n, dd))
    mses.append(float(((x - xh) ** 2).mean()))
    sq = ((x - xh) ** 2).sum(1)
    gll.append(float((-0.5 * dd * (LOG2PI + np.log(s2)) - sq / (2 * s2)).mean()))
    lll.append(float((-dd * np.log(2 * b) - np.abs(x - xh).sum(1) / b).mean()))
mses, gll, lll = np.array(mses), np.array(gll), np.array(lll)


def resid(xv, yv):
    A = np.vstack([xv, np.ones_like(xv)]).T
    c, *_ = np.linalg.lstsq(A, yv, rcond=None)
    return A @ c - yv


fig = make_subplots(
    rows=1, cols=2, horizontal_spacing=0.11,
    subplot_titles=(
        "Leg B — ℒ_div tracks Total Correlation (r = 0.9995)",
        "Leg A — MSE ⇔ log-likelihood: exact for Gaussian, broken for Laplace"))

fig.add_scatter(x=tcs, y=divs, mode="markers", name="ℒ_div vs TC (Gaussian latents)",
                marker=dict(size=8, color="#2F6F8F"), row=1, col=1)
A = np.vstack([tcs, np.ones_like(tcs)]).T
c, *_ = np.linalg.lstsq(A, divs, rcond=None)
fig.add_scatter(x=tcs, y=A @ c, mode="lines", name=f"best affine fit (slope {c[0]:.3f})",
                line=dict(color="#B07C2B", dash="dash"), row=1, col=1)
fig.update_xaxes(title_text="true Gaussian TC(F′) [nats]", row=1, col=1)
fig.update_yaxes(title_text="ℒ_div  (Eq. 3)", row=1, col=1)

fig.add_scatter(x=mses, y=resid(mses, gll), mode="markers+lines",
                name="Gaussian decoder — residual from affine",
                marker=dict(size=7, color="#2F6F8F"), row=1, col=2)
fig.add_scatter(x=mses, y=resid(mses, lll), mode="markers+lines",
                name="Laplace decoder (CONTROL) — residual",
                marker=dict(size=7, color="#A33B3B"), row=1, col=2)
fig.update_xaxes(title_text="reconstruction MSE (ℒ_rec)", row=1, col=2)
fig.update_yaxes(title_text="residual from best affine fit [nats]", row=1, col=2)

fig.update_layout(template="plotly_white", height=440,
                  legend=dict(orientation="h", y=-0.22),
                  title="Theorem 3.1 numerical audit — both legs, each with its control")
fig.write_html("results/claim2_theorem31.html", include_plotlyjs="cdn")

with open("results/claim2_theorem31.csv", "w", newline="") as fh:
    w = csv.writer(fh)
    w.writerow(["panel", "x", "y", "series"])
    for a, bb in zip(tcs, divs):
        w.writerow(["B_ldiv_vs_tc", f"{a:.6f}", f"{bb:.6f}", "ldiv"])
    rg, rl = resid(mses, gll), resid(mses, lll)
    for m, g, l in zip(mses, rg, rl):
        w.writerow(["A_residual", f"{m:.6f}", f"{g:.6e}", "gaussian"])
        w.writerow(["A_residual", f"{m:.6f}", f"{l:.6e}", "laplace_control"])

print(json.dumps({
    "ldiv_vs_tc_slope": float(c[0]),
    "ldiv_vs_tc_max_resid_nats": float(np.abs(A @ c - divs).max()),
    "gaussian_max_resid_nats": float(np.abs(resid(mses, gll)).max()),
    "laplace_control_max_resid_nats": float(np.abs(resid(mses, lll)).max()),
}, indent=2))
print("wrote results/claim2_theorem31.html + .csv")