1import json
  2import numpy as np
  3
  4np.set_printoptions(precision=7, suppress=True)
  5RNG = np.random.default_rng(7)
  6
  7
  8def random_hermitian(rng, d):
  9    z = rng.normal(size=(d, d)) + 1j * rng.normal(size=(d, d))
 10    h = (z + z.conj().T) / 2
 11    h -= np.trace(h).real * np.eye(d) / d
 12    return h / np.linalg.norm(h, ord=2)
 13
 14
 15def hermitian_function_and_derivatives(theta, ops, phi=np.tanh):
 16    """Matrix functional calculus and its exact Frechet derivatives."""
 17    H = sum(t * h for t, h in zip(theta, ops))
 18    lam, U = np.linalg.eigh(H)
 19    A = (U * phi(lam)) @ U.conj().T
 20    G = np.empty((len(lam), len(lam)))
 21    for a in range(len(lam)):
 22        for b in range(len(lam)):
 23            if abs(lam[a] - lam[b]) < 1e-9:
 24                G[a, b] = 1.0 - np.tanh(lam[a]) ** 2
 25            else:
 26                G[a, b] = (phi(lam[a]) - phi(lam[b])) / (lam[a] - lam[b])
 27    dAs = []
 28    for h in ops:
 29        hb = U.conj().T @ h @ U
 30        dAs.append(U @ (G * hb) @ U.conj().T)
 31    return A, dAs
 32
 33
 34def predictions_and_grad(theta, ops, rhos, y):
 35    A, dAs = hermitian_function_and_derivatives(theta, ops)
 36    yh = np.real(np.einsum('bij,ji->b', rhos, A))
 37    err = yh - y
 38    loss = float(np.mean(err * err))
 39    grad = np.array([2 * np.mean(err * np.real(np.einsum('bij,ji->b', rhos, dA))) for dA in dAs])
 40    return loss, yh, grad
 41
 42
 43def adam_fit(ops, rhos, y, steps=700, lr=0.06, seed=1):
 44    rng = np.random.default_rng(seed)
 45    theta = rng.normal(0, .12, len(ops))
 46    m = np.zeros_like(theta); v = np.zeros_like(theta)
 47    for t in range(1, steps + 1):
 48        _, _, g = predictions_and_grad(theta, ops, rhos, y)
 49        m = .9*m + .1*g; v = .999*v + .001*g*g
 50        theta -= lr * (m/(1-.9**t)) / (np.sqrt(v/(1-.999**t)) + 1e-8)
 51    return theta, predictions_and_grad(theta, ops, rhos, y)[0]
 52
 53
 54def classical_fit(features, y, steps=700, lr=.06, seed=1):
 55    rng = np.random.default_rng(seed)
 56    w = rng.normal(0, .12, features.shape[1])
 57    m = np.zeros_like(w); v = np.zeros_like(w)
 58    for t in range(1, steps + 1):
 59        e = features @ w-y; g = 2*features.T@e/len(y)
 60        m=.9*m+.1*g; v=.999*v+.001*g*g
 61        w -= lr*(m/(1-.9**t))/(np.sqrt(v/(1-.999**t))+1e-8)
 62    return w, float(np.mean((features@w-y)**2))
 63
 64
 65def main():
 66    d, ntrain, ntest = 4, 180, 180
 67    # Generic fixed Hermitian interactions avoid the special two-Pauli identity.
 68    ops = [random_hermitian(RNG, d) for _ in range(3)]
 69    # Commuting control: same parameter count, diagonal interactions.
 70    cops = [np.diag(RNG.uniform(-1, 1, d)).astype(complex) for _ in range(3)]
 71    target_theta = np.array([.85, -.65, .55])
 72
 73    def states(n):
 74        q = RNG.normal(size=(n,d)) + 1j*RNG.normal(size=(n,d))
 75        q /= np.linalg.norm(q, axis=1, keepdims=True)
 76        return q[:,:,None] * q[:,None,:].conj()
 77    rtrain, rtest = states(ntrain), states(ntest)
 78    At, _ = hermitian_function_and_derivatives(target_theta, ops)
 79    ytrain=np.real(np.einsum('bij,ji->b',rtrain,At))
 80    ytest=np.real(np.einsum('bij,ji->b',rtest,At))
 81
 82    # Core claim: divided-difference Frechet gradient agrees with finite differences.
 83    theta=np.array([.31,-.22,.17]); _,_,g=predictions_and_grad(theta, ops, rtrain[:40], ytrain[:40])
 84    eps=1e-5; fd=[]
 85    for j in range(3):
 86        tp=theta.copy(); tm=theta.copy(); tp[j]+=eps; tm[j]-=eps
 87        lp=predictions_and_grad(tp,ops,rtrain[:40],ytrain[:40])[0]
 88        lm=predictions_and_grad(tm,ops,rtrain[:40],ytrain[:40])[0]
 89        fd.append((lp-lm)/(2*eps))
 90    fd=np.array(fd)
 91    grad_rel_error=float(np.linalg.norm(g-fd)/(np.linalg.norm(fd)+1e-12))
 92    comms=[np.linalg.norm(ops[i]@ops[j]-ops[j]@ops[i]) for i in range(3) for j in range(i)]
 93    ccomms=[np.linalg.norm(cops[i]@cops[j]-cops[j]@cops[i]) for i in range(3) for j in range(i)]
 94
 95    results={}
 96    for name, model_ops in [('spectral_noncommuting',ops),('spectral_commuting',cops)]:
 97        theta_fit, train_loss=adam_fit(model_ops,rtrain,ytrain,seed=3)
 98        test_loss=predictions_and_grad(theta_fit,model_ops,rtest,ytest)[0]
 99        results[name]={'train_mse':train_loss,'test_mse':test_loss,'theta':theta_fit.tolist()}
100    def feats(r, model_ops):
101        return np.column_stack([np.ones(len(r))]+[np.real(np.einsum('bij,ji->b',r,h)) for h in model_ops])
102    wc, lc=classical_fit(feats(rtrain,ops),ytrain,seed=3)
103    results['classical_linear_energy']={'train_mse':lc,'test_mse':float(np.mean((feats(rtest,ops)@wc-ytest)**2)),'weights':wc.tolist()}
104    out={'gradient_check':{'analytic':g.tolist(),'finite_difference':fd.tolist(),'relative_error':grad_rel_error},
105         'commutator_frobenius':{'noncommuting':comms,'commuting':ccomms},'results':results,
106         'setup':{'seed':7,'train':ntrain,'test':ntest,'dimension':d,'parameters':3,'activation':'tanh','target':'tanh(sum(theta_j H_j))'}}
107    print(json.dumps(out, indent=2))
108
109if __name__=='__main__': main()