Affine directional Sobolev regularizer

Integrate the method at the loss level for a classifier or regressor whose input is a vectorized image, tabular feature vector, or embedding. On each selected regularization step, draw K Gaussian vectors r_k in R^n and normalize them as \(\xi_k=r_k/(\|r_k\|_2+10^{-12})\). Draw a normalized output projection u, or apply the penalty independently to selected output coordinates. For every minibatch input x, compute the scalar output \(s(x)=u^\top f_\theta(x)\), then obtain the directional derivative \(d_{x,k}=\xi_k^\top\nabla_xs(x)\) with reverse-mode autograd or a Jacobian-vector product. Compute \(a_k=(|B|^{-1}\sum_x|d_{x,k}|^p+\varepsilon)^{1/p}\), combine them as \(E=(K^{-1}\sum_ka_k^{-n})^{-1/n}\), and backpropagate \(L_{task}+\lambda E\). Start with p=2 or p=4, K=4, and epsilon=1e-8; clip or log-transform E if occasional nearly zero directional norms create exploding gradients. Also test a scale-normalized version \(E/(|B|^{-1}\sum_x\|f_\theta(x)\|_2^2+\varepsilon)^{1/2}\). The paper provides the spherical inverse-power aggregation; minibatch integration, Monte Carlo directions, epsilon stabilization, and output projection are approximations to make it usable in finite networks. First run on MNIST or CIFAR-10 with a 3-layer MLP or small CNN, comparing no Jacobian penalty, standard squared input-gradient regularization, and this penalty at matched model size and training budget. Measure clean accuracy, corruption accuracy, loss variance, and sensitivity after applying random anisotropic feature scaling or an invertible linear input transform. Success is equal or better accuracy and robustness with stable training and at most 1.5x wall-clock overhead.