Variable-Projection Adaptive Fourier Layer

Integrate this as the first layer of a coordinate MLP, PINN, or neural operator that maps coordinates x to a field value. Maintain M complex directions d_m and use features z_m(x)=exp(i k d_m^T x); concatenate real and imaginary parts with ordinary MLP features. For the MVP, make the Fourier-feature output linear: c is solved analytically at every outer step, while theta is optimized by Adam or L-BFGS. At each iteration, sample coordinates x_i and targets or PDE residual equations y_i; construct B_theta with B[i,m]=exp(1j*k*dot(d_m,x_i)); form A=B.conj().T@W@B+lambda*I and b=B.conj().T@W@y; solve A c=b using Cholesky rather than an explicit inverse; compute Phi; backpropagate through the linear solve to theta; and update theta. Every K steps compute G=B.conj().T@W@B, eigendecompose G, and prune features whose associated subspace has eigenvalue below tau*lambda_max(G). Optionally birth a new direction at the largest residual Fourier correlation. Parameterize real propagating directions in two dimensions by d=(cos(theta),sin(theta)); for evanescent waves allow theta=a+ib and evaluate complex sine and cosine in double precision. The paper supplies the residual minimization, variable projection, quadratic-growth intuition, and Gramian cutoff; estimate W, lambda, tau, and local curvature empirically. First test on a 2D Helmholtz field made from M=4, 8, or 16 plane waves, comparing a fixed random Fourier-feature PINN and a SIREN at equal feature count. Measure relative L2 field error, boundary/PDE residual, number of retained features, and wall-clock time. Success is lower error with fewer features or at least 2x faster residual reduction at equal parameter count; failure is inability to recover planted directions or unstable Gramian pruning.