Hierarchical B-spline sparse-grid front end

Integrate this as the input embedding immediately before the first MLP layer in a coordinate network, PINN, neural operator, or point-cloud model. Normalize each input \(x\) to \([0,1]^d\), choose spline degree \(p\), maximum level \(n\), and define \(\mathcal{I}_n=\{(\ell,i):\ell\in\mathbb{N}_0^d, |\ell|_1\le n, i\text{ is valid at level }\ell\}\). Construct compactly supported tensor-product B-splines and evaluate only active basis functions. The forward pass is features = [phi_k(x) for k in I_n]; z = Linear(features); output = MLP(z). For initialization from supervised samples \((x_s,t_s)\), form b[j] = sum_s w_s * phi_j(x_s), estimate M[j,k] = mean_x phi_j(x) * phi_k(x) using quadrature or uniform samples, and solve M c = b; use \(c\) to initialize the first linear layer or a coefficient field. Freeze the basis for the first training phase, then optionally unfreeze spline amplitudes while keeping knot locations fixed. Compute the sparse basis structure and Galerkin matrix from the paper's mathematics; estimate \(h\), effective sample count \(N\), and approximation error empirically. First run 2D and 3D coordinate regression on a smooth manufactured function and an anisotropic mixed-regularity function. Compare against dense tensor B-spline features and Fourier features at matched parameter counts and FLOPs. Sweep \(n\) and \(p\). Success is at least 2x fewer embedding parameters at equal validation error, or lower error at equal parameters, with stable optimization. Also test whether validation error follows the predicted tradeoff: increasing \(n\) initially lowers bias and eventually produces diminishing returns or noise amplification.