Certified nonconvex p-shrinkage

Integrate this into a proximal-gradient optimizer for a selected parameter tensor, such as all weights of a small MLP or only the router weights of an MoE layer. At iteration t, compute the ordinary gradient step z=W_t-eta_t*grad_W L(W_t), then apply the scalar proximal map independently to every entry of z with lambda=eta_t*alpha, where alpha is the l_p regularization coefficient. For each z_i, let r=abs(z_i), compute u_th=[2*lambda*(1-p)]**(1/(2-p)) and r_th=((2-p)/(2*(1-p)))*u_th; if r<r_th set W_new_i=0. Otherwise compute s=lambda*p*r**(p-2). If s is below a configurable series limit, initialize u=r*(1-s-sum_{n=2}^N c_n*s**n), evaluating c_n with log-Gamma functions, and accept it only when abs(h(u)-r)<=epsilon. Otherwise set lo=u_0=[lambda*p*(1-p)]**(1/(2-p)) and hi=r, then perform 30-50 bisection iterations on h(u)-r for the larger-root branch; use Newton steps only when they stay inside the bracket. The mathematical quantities are the threshold, stationary-root equation, and series coefficients. Estimate truncation reliability by comparing orders N and 2N; certify root accuracy using residual and bracket width. Start with a 2-layer width-512 MLP on MNIST or CIFAR-10. Compare AdamW plus magnitude pruning, soft-threshold proximal l1, and direct subgradient l_p training at matched sparsity and FLOPs. Measure validation accuracy, exact zero fraction, proximal-call time, wall-clock time, and loss descent. Success means higher accuracy at the same zero fraction, or equal accuracy with fewer pruning operations, with no threshold oscillations and residuals below 1e-6.