Differentiable Hankel PSD regularizer

(1) Integration point: add the regularizer to a conditional VAE, scalar density network, or world model whose output is a distribution over a bounded scalar. Start with a beta head: the network emits p_0 and q_0, then compute p=softplus(p_0)+epsilon and q=softplus(q_0)+epsilon. A general version emits moments or free cumulants a_0 through a_{2N_max-2}. (2) Pseudocode: compute s=p+q, r=p*q, and g=2*s**3*(s+1)+r*(s**3-7*s**2-16*s-12). Add w_g*softplus(-g/tau_g) to the task loss. For every selected N, construct H[i,j]=a[i+j], symmetrize it with H=(H+H.T)/2, compute its smallest eigenvalue, and add w_N*softplus((epsilon-lambda_min)/tau). Backpropagate through the eigendecomposition and the network. (3) The exact quantities from the paper are the beta certificate g and the Hankel PSD principle. Estimate the moments or cumulants from the model output; tune epsilon, temperatures, weights, and N_max empirically. If eigenvalue gradients are unstable, use a jittered Cholesky failure penalty or random-vector violations max(0,-v.T@H@v). (4) First experiment: train a small conditional VAE on synthetic mixtures of beta distributions. Compare likelihood-only training, a beta head without constraints, and the constrained head using N=2,3,4. Measure held-out likelihood, moment calibration, negative-eigenvalue frequency, and gradient stability. Success is a lower invalid-distribution rate and better high-order moment calibration at similar likelihood and compute.