Operator-spline theory: consolidated research manuscript

Rendered archival TeX

This is an HTML reading rendition of the local TeX record. Mathematical notation is rendered with KaTeX; archived figures are included when their source assets are part of this collection.

documentclass[11pt]article

usepackage[margin=1in]geometry usepackage(amsmath,amssymb,amsthm,mathtools) usepackage(bm) usepackage(booktabs) usepackage(graphicx) usepackage(longtable) usepackage(hyperref) usepackage(enumitem) usepackage(xcolor) usepackage(float)

hypersetup( colorlinks=true, linkcolor=blue!50!black, citecolor=blue!50!black, urlcolor=blue!50!black )

newtheorem(definition)Definition newtheorem(proposition)Proposition newtheorem(remark)Remark

newcommand(\R)mathbb(R) newcommand(\Z)mathbb(Z) newcommand(\C)mathbb(C) newcommand(\dd)mathrm(d) newcommand(\calL)mathcal(L) newcommand(\calN)mathcal(N) newcommand(\bs)boldsymbol newcommand(\A)mathbf(A) newcommand(\G)mathbf(G) newcommand(\I)mathbf(I) newcommand(\cvec)mathbf(c) newcommand(\xvec)mathbf(x) newcommand(\yvec)mathbf(y) newcommand(\wvec)mathbf(w) newcommand(\zvec)mathbf(z) newcommand(\uvec)mathbf(u) DeclareMathOperator*argminarg min

title(Operator-Spline Neural Representations: \\ Continuous-Domain Operator Bases for Efficient Physical Learning) author(OSNR Project Notes) date(\today)

begin(document) maketitle

tableofcontents

Introduction

Coordinate networks have become a standard tool for representing continuous signals and fields. A typical implicit neural representation (INR) maps coordinates to field values through a multilayer perceptron, $$f_\theta : \xvec \mapsto y,$$ where the weights $\theta$ are optimized by stochastic gradient descent. SIRENs improve high-frequency representation by using sinusoidal activations, while PINNs add differential-equation residuals to the training objective.

The OSNR thesis is that this approach is structurally misaligned for physical systems governed by known operators. It is inspired by the operator-based spline signal-processing program of Unser, Blu, Vetterli, and collaborators [unser1993bspline1,unser1993bspline2,unser2005cardinal1,unser2005cardinal2,vetterli2002fri]. If a field is constrained by $$L\{s\} = r,$$ then the representation should be built from the operator $L$ itself. The role of learning or numerical inversion should be reduced to coefficient recovery in an already appropriate continuous function space.

OSNR therefore replaces black-box nonlinear layers by operator-matched spline dictionaries. The spline coefficients play the role of the latent representation. The architecture is not a generic neural network with a different activation; it is a continuous-domain inverse-problem engine presented in neural-representation form, aligned with continuous-domain inverse-problem representer theorems [gupta2018continuous,debarre2019hybrid,debarre2021composite]. The same principle also suggests a deterministic alternative to the trunk side of branch-trunk neural operator models such as DeepONet and MIONet [lu2021deeponet,jin2022mionet]: the coordinate-to-field map can be a spline synthesis operator with exact calculus rather than a learned MLP.

The empirical record tests this thesis rather than establishing universal superiority. S[sec:grown-topology] covers closed-form identification, grown-topology control and continual-learning comparisons; S[sec:matched-vs-sindy] compares operator-matched identification with SINDy and neural ODEs, including a negative blind-discovery control. S[sec:pde-vs-fno] studies PDE comparisons with neural operators, regime changes and partially known physics. Their measured gains depend on the information supplied, comparator, accuracy and computational accounting: knowing an operator family is not by itself a sufficient condition for an advantage. The later exact battery-response experiment, for example, admits a nearly equally fast training-free classical control; the real battery-aging fit loses to its matched polynomial comparator. S[sec:ssp-view] develops the sparse-stochastic-process connection under its stated assumptions, not a universal statistical dominance theorem. S[sec:rsi] records recursive-improvement experiments, while subsequent claim audits distinguish pooled fitting statistics, protected function supports, immutable archives and retained task performance. None alone establishes unrestricted self-improvement without forgetting or drift. The constructive result is a tested representation and calculus toolbox; independent-task benefit against strong controls remains the capability test.

Mathematical Background

Cardinal spline representation

The classical cardinal spline model represents a continuous signal as [unser1993bspline1,unser1993bspline2] $$s(t) = \sum_{k \in \Z} c[k] \beta(t-k),$$ where $\beta$ is a spline generator and $c[k]$ are discrete coefficients. The essential point is that a continuous function space is controlled by a discrete sequence. This makes exact digital processing of continuous objects possible, provided the generator is stable and the correct coefficient-domain filters are used.

Linear differential operators and null spaces

Let $L$ be a constant-coefficient differential operator with characteristic roots or poles $$\bs{\alpha} = (\alpha_1,\ldots,\alpha_N).$$ The null space is $$\calN_L = \{s : L\{s\}=0\} = \mathrm{span}\{e^{\alpha_m t}\}_{m=1}^N$$ with polynomial factors included for repeated roots. For the Helmholtz operator $$L = \frac{\dd^2}{\dd x^2} + k^2,$$ the poles are $\alpha=\pm j k$, and the null space consists of sinusoidal modes.

Cardinal exponential splines

Cardinal exponential splines are compactly supported spline generators matched to $\bs{\alpha}$ [unser2005cardinal1]. Their Fourier-domain form is $$\widehat{\beta}_{\bs{\alpha}}(\omega) = \prod_{m=1}^N \frac{1 - e^{\alpha_m - j\omega}}{j\omega-\alpha_m}.$$ Integer shifts of $\beta_{\bs{\alpha}}$ generate a stable spline space under appropriate Riesz conditions. The basis reproduces exponential polynomials and supports exact operator calculus in the coefficient domain.

Exponential Hermite splines

Cardinal E-splines attach the operator to a scalar coefficient stream. For higher-order boundary value problems, OSNR requires a vector-valued cardinal system whose coefficients store not only function values but also derivatives. The second-order exponential-polynomial Hermite generator of Schmitter, Badoual, Uhlmann, Fageot, and Unser [schmitter2016hermite] provides exactly this structure.

Let $$\Phi(t)= \begin{bmatrix} \phi_0(t) & \phi_1(t) & \phi_2(t) \end{bmatrix}^{\top}.$$ The Hermite spline expansion is $$f(t) = \sum_{k\in\Z} \left( c_0[k]\phi_0(t-k) +c_1[k]\phi_1(t-k) +c_2[k]\phi_2(t-k) \right),$$ with coefficient vectors $$\mathbf{c}[k]= \begin{bmatrix} c_0[k] \\ c_1[k] \\ c_2[k] \end{bmatrix} = \begin{bmatrix} f(k) \\ f'(k) \\ f''(k) \end{bmatrix}$$ for exactly interpolated data. The interpolation conditions are $$\phi_p^{(r)}(k)=\delta_{pr}\delta_k, \qquad p,r\in\{0,1,2\},\quad k\in\Z.$$ On the interval $[0,1]$, each channel has the exponential-polynomial form $$\phi_p(t) = A_p+B_p t+C_p t^2+D_p t^3 +E_p e^{j\omega_0 t}+F_p e^{-j\omega_0 t}, \qquad \omega_0=\frac{2\pi}{M}.$$ The constants $A_p,\ldots,F_p$ are determined by the six Hermite endpoint constraints $$\phi_p^{(r)}(0)=\delta_{pr}, \qquad \phi_p^{(r)}(1)=0, \qquad r=0,1,2.$$ The negative branch is fixed by the Hermite parity $$\phi_0(-t)=\phi_0(t),\qquad \phi_1(-t)=-\phi_1(t),\qquad \phi_2(-t)=\phi_2(t),$$ so the generator is compactly supported on $[-1,1]$. This construction is $C^2$, reproduces polynomials up to cubic degree, and reproduces the trigonometric modes $\sin(\omega_0 t)$ and $\cos(\omega_0 t)$ through the coefficient samples of the function and its first two derivatives.

For OSNR Tier 3, the crucial change is semantic: Dirichlet, Neumann, and curvature boundary data become direct coefficient assignments. A clamped boundary at knot $k_b$, for example, is imposed by setting $$\mathbf{c}[k_b] = \begin{bmatrix} 0 \\ 0 \\ 0 \end{bmatrix},$$ rather than by adding a soft loss term. The same vector-valued structure also produces the block-circulant Hermite Gram system derived in Appendix [app:hermite-block-gram].

OSNR Architecture

Master core apparatus

The current OSNR research code separates into three production tiers. Each tier has a different admissible function space, solver topology, and verification target.

Tier 1: deterministic operator-bound systems

Tier 1 targets smooth physical fields that lie in, or close to, the null space of a known linear operator. Examples include Helmholtz waves, damped oscillators, and linear constant-coefficient PDE components. The solver pipeline is:

1. identify $L$ and its poles $\bs{\alpha}$;

2. construct calibrated E-spline or null-space dictionaries;

3. recover coefficients through a stable linear solve or circulant FFT inversion;

4. compute differential quantities through operator identities rather than autograd.

Tier 2: adaptive sparse continua

Tier 2 targets composite fields with sparse discontinuities or shocks: $$s = s_{\mathrm{smooth}} + s_{\mathrm{sparse}}.$$ Uniform grids are not sufficient for non-bandlimited discontinuities at sub-grid coordinates. The sparse tier must first identify finite-rate innovations, then adapt the dictionary to those coordinates:

1. estimate innovation locations using FRI or matrix-pencil methods;

2. snap sparse knots to the recovered coordinates;

3. solve a cross-Gram-coupled sparse-plus-smooth inverse problem;

4. debias with scale-invariant ridge stabilization.

Tier 3: higher-order neural operators

Tier 3 targets families of PDE solutions rather than a single fitted field. Neural operators such as DeepONet and MIONet learn maps between function spaces by pairing branch networks, which encode input functions or boundary data, with trunk networks, which encode query coordinates [lu2021deeponet,jin2022mionet]. Spline-PINN shows a related but distinct path: a CNN predicts Hermite spline coefficients, and a continuous Hermite spline layer evaluates PDE residuals without finite-difference losses [wandel2022splinepinn]. OSNR adopts the continuous Hermite idea but removes the black-box coordinate trunk where the governing operator and boundary calculus are known.

Let $$\Phi(t) = \begin{bmatrix} \phi_0(t) & \phi_1(t) & \phi_2(t) \end{bmatrix}^{\top}$$ be a second-order Hermite generator compactly supported on $[-1,1]$. The generalized higher-order Hermite construction of Schmitter, Badoual, Uhlmann, Fageot, and Unser stores value, slope, and curvature data at each knot and interpolates all three channels exactly [schmitter2016hermite]: $$\phi_0^{(r)}(k)=\delta_{r0}\delta_k,\qquad \phi_1^{(r)}(k)=\delta_{r1}\delta_k,\qquad \phi_2^{(r)}(k)=\delta_{r2}\delta_k, \qquad r=0,1,2.$$ The synthesized field is $$f(t) = \sum_{k\in\Z} c_0[k]\phi_0(t-k) + c_1[k]\phi_1(t-k) + c_2[k]\phi_2(t-k),$$ where $$\mathbf{c}[k] = \begin{bmatrix} c_0[k] \\ c_1[k] \\ c_2[k] \end{bmatrix} = \begin{bmatrix} f(k) \\ f'(k) \\ f''(k) \end{bmatrix}$$ for exactly interpolated data. The Hermite generator in [schmitter2016hermite] is piecewise polynomial-exponential, $C^2$, compactly supported, and reproduces cubic polynomials as well as trigonometric functions. This is the missing boundary mechanism for a neural-operator OSNR tier: Dirichlet, Neumann, and curvature constraints are coefficient assignments, not soft loss penalties.

For a PDE solution operator $$\mathcal{S}: u \mapsto v,$$ the branch-side network or analytic encoder should output Hermite coefficient tensors $$\{\mathbf{c}[k;u]\}_{k\in\Z},$$ while the trunk side is replaced by deterministic Hermite synthesis. In multidimensional domains, tensor products of one-dimensional Hermite generators give mixed value/derivative channels, matching the construction used by Spline-PINN for continuous PDE residuals [wandel2022splinepinn]. Unlike Spline-PINN, OSNR evaluates continuous energies and cross-correlations in coefficient space through the block-Gram calculus derived in Appendix [app:hermite-block-gram], avoiding Monte Carlo residual quadrature whenever the operator and boundary model admit exact inner products.

Continuous-Discrete Calibration

The most important implementation condition is the cardinal coordinate map, which preserves the shift-invariant structure required by spline filtering calculus [unser2005cardinal1,unser2005cardinal2]. $$v_k(x) = \frac{x}{T} - k.$$ Here $T$ is the physical knot spacing. For a domain $[0,D]$ and a compact generator of support order $q$, a stable finite dictionary uses $$T = \frac{D}{M-q}.$$ Then $$v_k(x+T) = v_k(x)+1,$$ so one physical knot step maps to one cardinal interval.

This explains why a derivative residual can be zero while reconstruction fails. If the derivative dictionary is defined algebraically by $\A_{d2}=-k^2 \A$, then the PDE residual $$\A_{d2}\cvec + k^2 \A \cvec$$ is identically zero regardless of whether $\A$ is a well-conditioned reconstruction basis.

2D tensor-product Hermite expansion

The 2D Tier 3 engine is obtained by tensorizing the one-dimensional second-order Hermite streams. Let $$h_i^x(x),\qquad h_j^y(y),\qquad i,j\in\{0,1,2\},$$ denote the value, slope, and curvature Hermite generators along the two axes. The tensor-product basis functions are $$h_{i,j}(x,y) = h_i^x(x)h_j^y(y), \qquad i,j\in\{0,1,2\}.$$ For a grid node $(k,\ell)$, OSNR stores a nine-stream local state $$\mathbf{c}_{k,\ell} = \begin{bmatrix} f & \partial_x f & \partial_{xx}f & \partial_y f & \partial_{xy}f & \partial_{xxy}f & \partial_{yy}f & \partial_{xyy}f & \partial_{xxyy}f \end{bmatrix}_{(k,\ell)}^{\top}.$$ The synthesized field is $$f(x,y) = \sum_{k,\ell} \sum_{i=0}^{2}\sum_{j=0}^{2} c_{i,j}[k,\ell]\, h_i^x\!\left(\frac{x}{T_x}-k\right) h_j^y\!\left(\frac{y}{T_y}-\ell\right).$$ This formula is the tensor-product analogue of the Schmitter et al. Hermite generator and the multidimensional counterpart of the periodic inner-product calculus of Badoual, Schmitter, and Unser.

The practical consequence is that the spatial calculus of 2D physical fields is a forward coefficient operation. In the stream-function formulation for incompressible flow, $$v_x=\partial_y a_z,\qquad v_y=-\partial_x a_z,$$ and therefore $$\nabla\cdot \mathbf{v} = \partial_x\partial_y a_z-\partial_y\partial_x a_z = 0$$ up to the commutation error of the discrete finite-difference stencil. The nonlinear transport term is evaluated as $$(\mathbf{v}\cdot\nabla)\mathbf{v} = \begin{bmatrix} v_x\partial_x v_x+v_y\partial_y v_x\\ v_x\partial_x v_y+v_y\partial_y v_y \end{bmatrix},$$ and viscous diffusion as $$\Delta\mathbf{v} = \begin{bmatrix} \partial_{xx}v_x+\partial_{yy}v_x\\ \partial_{xx}v_y+\partial_{yy}v_y \end{bmatrix}.$$ All operators appearing in these expressions are evaluated by shift-invariant finite-difference ladders on the Hermite coefficient streams during the forward pass. No backward-mode automatic differentiation tape is constructed; the measured PyTorch autograd graph allocation in all Tier 3 validations is therefore $0.00$ bytes.

Autograd-Free Differential Calculus

For an exponential spline with pole vector $\bs{\alpha}$, applying a first-order operator $(D-\alpha_m)$ reduces the order of the spline [unser2005cardinal1,delgadogonzalo2012exponential]: $$(D-\alpha_m)\beta_{\bs{\alpha}}(t) = \beta_{\bs{\alpha}\setminus \alpha_m}(t) - e^{\alpha_m} \beta_{\bs{\alpha}\setminus \alpha_m}(t-1).$$ This is the finite-difference ladder. Differential fields can be evaluated by filtering coefficients or by applying deterministic lower-order dictionary maps. For the Helmholtz null space, $$\frac{\dd^2}{\dd x^2} s(x) = -k^2 s(x)$$ inside the smooth spans, with boundary innovations handled separately.

This removes the need for backward-mode automatic differentiation in Tier 1 PDE residual evaluation. The differential operator is encoded in the basis and coefficient algebra.

Circulant and FFT Solvers

When the dictionary is shift-invariant and periodized, the Gram matrix is circulant. This is the coefficient-domain form of the exact periodic inner-product calculus for spline curves and functions [badoual2016inner,badoual2018periodic]: $$\G = \begin{bmatrix} g_0 & g_{M-1} & \cdots & g_1 \\ g_1 & g_0 & \cdots & g_2 \\ \vdots & \vdots & \ddots & \vdots \\ g_{M-1} & g_{M-2} & \cdots & g_0 \end{bmatrix}.$$ Such matrices are diagonalized by the discrete Fourier transform: $$\G = \mathbf{F}^{-1} \Lambda \mathbf{F}.$$ Solving $\G \cvec = \mathbf{b}$ reduces to $$\widehat{\cvec}[\ell] = \frac{\widehat{\mathbf{b}}[\ell]}{\lambda_\ell}.$$ This converts $O(M^3)$ dense solves into $O(M\log M)$ FFT operations.

Tomographic Radiance Fields as Spline Inverse Problems

The failed direct ray-kernel DL3DV experiment clarifies an important modeling boundary. A nontrivial view-synthesis scene is not naturally a smooth map from ray origin and direction to RGB. The physically shared quantity is a latent field in space, observed through line or ray measurements. The spline tomography literature gives the appropriate replacement model: represent the unknown continuous field by shifted basis functions, push those basis functions through the forward projector, and solve the coefficient inverse problem with fast adjoint and normal operators [nilchian2013fast,mccann2016fast,donati2018multiscale,haouchat2025generalized]. Jin et al. make the complementary point that inverse problems whose normal operators are convolutional admit physics-aware direct inversions before any learned artifact-removal stage [jin2017deep]; in OSNR, the direct inverse is the primary object, and any neural residual must remain secondary.

For a first linearized density or opacity stage, write $$\sigma(\mathbf{x})=\sum_{\mathbf{k}} a[\mathbf{k}]\,\varphi_\sigma(\mathbf{x}-\Lambda\mathbf{k}), \qquad \mathbf{g}=H\mathbf{a}+\boldsymbol{\eta},$$ where $H$ samples line or ray integrals of the spline density field. The regularized inverse update is $$\mathbf{a}^{\star} = \arg\min_{\mathbf{a}} \frac12\|W^{1/2}(H\mathbf{a}-\mathbf{g})\|_2^2 +\lambda R(\mathbf{a}),$$ with $W$ encoding ray confidence or frequency reliability, as in weighted phase-retrieval formulations [bostan2016variational]. For quadratic $R(\mathbf{a})=\|L\mathbf{a}\|_2^2$, the normal equation is $$(H^\top W H+\lambda L^\top L)\mathbf{a} = H^\top W\mathbf{g}.$$ The computational opportunity is the McCann–Donati normal-operator identity. For shift-invariant basis functions and locally stationary projection blocks, $H^\top H$ acts as a discrete convolution: $$(H^\top H\mathbf{a})[\mathbf{k}] = (\mathbf{a}*\mathbf{r})[\mathbf{k}],$$ where $\mathbf{r}$ is the sampled autocorrelation of the projected basis function. Thus the expensive repeated normal-operator application inside conjugate gradients or ADMM becomes a Fourier-domain multiplication. Multiscale basis functions then provide a controlled coarse-to-fine path that is robust to pose and angular uncertainty [donati2018multiscale].

Fast forward projection and sparse acquisition variants can be imported from the same lineage: Arcadu et al. use Fourier regridding with minimal oversampling for efficient forward projectors, while Donati et al. show how randomized STEM sampling can be coupled to regularized tomographic recovery [arcadu2016forward,donati2017compressed]. This does not make full NeRF rendering linear. The volume-rendering equation contains transmittance $$C(r)=\int T(t)\sigma(r(t))c(r(t),\mathbf{d})\,dt, \qquad T(t)=\exp\!\left(-\int_0^t\sigma(r(s))\,ds\right),$$ so exact RGB fitting remains nonlinear in $\sigma$. The OSNR route is therefore staged: first recover coarse density/support with a tomographic spline inverse solve; then refine the density multiscale; then solve color/radiance coefficients on the recovered support; and finally apply visibility-weighted nonlinear corrections. Positivity, support constraints, total-variation, Hessian, and sparse-innovation priors enter naturally through the constrained ADMM machinery developed for spline tomography [nilchian2013constrained,nilchian2015spline].

The controlled runner apps\_industrial\_breakthrough/spline\_tomographic\_radiance\_solver.py validates only the linearized operator claim. It builds a periodic synthetic density field, samples $48$ discrete projection directions, constructs $H^\top H$ explicitly once from a delta impulse, and then replaces all subsequent normal-operator applications by FFT convolution. On a $96\times96$ field, the convolutional normal operator matches explicit $H^\top H$ with relative error $2.4462\times10^{-7}$. One explicit normal-operator application costs $2.4956$ ms, while the FFT version costs $0.0969$ ms. Solving the same ridge-regularized inverse problem by conjugate gradients takes $188.1746$ ms with explicit normals and $4.9666$ ms with FFT normals, yielding a reconstruction PSNR of $22.3578$ dB under noisy sparse projections. This is not yet a NeRF result; it is the isolated mathematical validation that the spline-tomographic normal operator can be diagonalized as the literature predicts.

The next controlled runner, apps\_industrial\_breakthrough/spline\_ray\_operator\_validation.py, validates the more fundamental Haouchat-style requirement: the forward ray operator and its adjoint must be matched before any real-scene radiance experiment is meaningful. On a $36\times36$ coefficient grid with $2688$ parallel rays, the script constructs two explicit small operators for auditability: a pixel basis and a quadratic tensor-product spline basis. The spline data are generated by the spline operator itself, and both models solve the same noisy inverse problem with conjugate gradients. The adjoint identity $\langle Hc,p\rangle=\langle c,H^\top p\rangle$ holds to $1.5672\times10^{-15}$ relative error for the spline operator and $2.7427\times10^{-16}$ for the pixel operator. At the same coefficient count, the spline inverse reconstructs the continuous rendered target at $54.0869$ dB, while the pixel model reaches only $27.0125$ dB. This is an intentionally controlled operator test: it proves that the coefficient-domain ray basis and adjoint are now correctly formulated, not that the full DL3DV visibility problem is solved.

The basis choice itself is not incidental. Following the exponential-spline construction of Delgado-Gonzalo, Thevenaz, and Unser [delgadogonzalo2012exponential], a one-dimensional cardinal exponential B-spline associated with poles $\boldsymbol{\alpha}=(\alpha_1,\ldots,\alpha_N)$ has Fourier-domain form $$\widehat{\beta}_{\boldsymbol{\alpha}}(\omega) = \prod_{m=1}^{N} \frac{1-\exp(\alpha_m-j\omega)}{j\omega-\alpha_m}.$$ The two-dimensional smooth tier then uses the tensor-product generator $$\varphi_{\boldsymbol{\alpha}_x,\boldsymbol{\alpha}_y}(x,y) = \beta_{\boldsymbol{\alpha}_x}(x)\, \beta_{\boldsymbol{\alpha}_y}(y),$$ so the basis can reproduce the local modes implied by the operator rather than merely interpolate samples. The controlled sweep apps\_industrial\_breakthrough/exponential\_spline\_basis\_sweep.py tests this precision-first hypothesis on the same $1920$ rays and $784$ coefficients while changing only the tensor-product basis. The data are generated by a damped-harmonic exponential spline with poles $(-\lambda,-\lambda+j\omega,-\lambda-j\omega)$, and all candidate operators reuse the identical ray geometry. The matched damped-harmonic basis reaches $56.3837$ dB, compared with $52.1420$ dB for a monotone exponential-decay basis, $50.0620$ dB for an undamped harmonic basis, $49.9084$ dB for a quadratic polynomial spline, and $28.7579$ dB for a pixel box basis. This isolates the important design rule for the next radiance-field stage: the highest precision should come from tensor-product operator splines whose poles are matched to the expected local dynamics, with sparsity and compression applied after the operator basis is correct.

The follow-up runner apps\_industrial\_breakthrough/exponential\_spline\_pole\_sweep.py turns this from a hand-picked basis comparison into a deterministic pole-selection problem. It fixes the target field, ray geometry, coefficient count, noise level, ridge parameter, and conjugate-gradient budget, then sweeps damped-harmonic exponential splines of orders $2$, $3$, and $4$ over $\lambda\in\{0.18,0.28,0.35,0.42,0.56,0.72\}$ and periods $\{8,10,12,16\}$. The target is generated by the order-$3$ pole set $(-0.42,-0.42+j2\pi/10,-0.42-j2\pi/10)$. The pole sweep correctly ranks that matched operator first at $55.2077$ dB on $1280$ rays and $576$ coefficients. The best order-$4$ candidate reaches $52.4859$ dB, and the best order-$2$ candidate reaches only $42.1900$ dB. This is the first automated evidence that pole selection is a meaningful OSNR model-selection axis: increasing support/order blindly does not dominate, while matching the operator poles controls reconstruction precision.

The support/regularity runner apps\_industrial\_breakthrough/exponential\_spline\_support\_regularization\_sweep.py then isolates the opposite regime: a field generated by the shortest first-order Green-matched spline with pole $\alpha=-0.42$. The one-pole basis has support length $1$ and no continuity guarantee, but it exactly matches the local Green mode. It reaches $65.7918$ dB with operator density $0.0362$ and an average of $20.86$ active coefficients per ray. The smooth order-$3$ mixed basis $[0,0,\alpha]$ reaches only $27.1299$ dB and has density $0.1084$ with $62.44$ active coefficients per ray. This confirms a second design rule: if the modeled object is a Green response or sparse innovation, the shortest matched spline can be both more accurate and more localized than a smoother high-order basis. Regularity should be introduced because the signal class requires it, not by default.

Finally, apps\_industrial\_breakthrough/exponential\_spline\_basis\_selection\_map.py evaluates the full pole-multiset selection problem. It tests seven target regimes against fourteen candidate bases: first-order Green response, repeated real poles, two distinct real poles, zero-augmented smooth operators, damped oscillators, pure polynomial splines, and support-$4$ mixed bases. Across all regimes, the exact pole multiset ranks first. This is the strongest evidence so far that OSNR basis design should be formulated as operator pole selection rather than degree selection. Order controls support and regularity; pole multiplicity and location control the reproduced null-space modes. Higher order improves asymptotic approximation power for smooth functions, but it is not a substitute for matching the operator that generated the signal.

The final synthetic step, apps\_industrial\_breakthrough/exponential\_spline\_operator\_inference.py, removes access to rendered target PSNR during model selection. Each candidate basis is fitted on $75\%$ of the rays and scored on held-out rays using a normalized validation residual plus small support and density penalties. This exposes a practical distinction between generative pole matching and predictive operator selection. Under sparse-ray training, the oracle PSNR basis is sometimes a smoother support-$4$ model rather than the exact generating basis, because the added regularity improves interpolation across unobserved rays. The validation score still selects the oracle basis in four of seven regimes and stays within $0.2039$ dB of oracle in all cases, with mean PSNR loss $0.0465$ dB. Thus the operational rule becomes: use the differential operator poles as the first prior, then choose among nearby pole augmentations by held-out measurement prediction rather than training residual.

We then stress-test the same idea in a mixed local-operator field with apps\_industrial\_breakthrough/exponential\_spline\_local\_operator\_adaptation.py. The synthetic field assigns different pole multisets to different spatial regions and compares a single global basis, a held-out-ray local block selector, and an oracle local label map. The oracle local solve reaches $33.5291$ dB, while the best global held-out model reaches $28.6559$ dB. This proves that local operator adaptation has substantial headroom. However, a blind ray-only greedy block selector reaches only $28.8074$ dB with $25\%$ block-label accuracy. Global line measurements make small independent block labels weakly identifiable unless the selection objective includes stronger spatial priors, localized measurements, or joint segmentation/coefficient optimization.

The sparse-innovation remedy is tested in apps\_industrial\_breakthrough/exponential\_spline\_fri\_region\_adaptation.py. Instead of selecting arbitrary blocks, the method first reconstructs a global proxy, computes derivative-energy profiles, extracts sparse FRI-style transition proposals, expands them into a small boundary lattice, and then jointly scores region geometry and region pole choices on held-out rays. The raw derivative peaks locate approximate boundaries at $x=(-1.4149,1.4149)$ and $y=3.0319$; held-out refinement moves them to $x=(-3.4149,3.4149)$ and $y=4.0319$, yielding $100\%$ region-label agreement at the coefficient-grid resolution. The resulting FRI-region adaptive model reaches $34.2000$ dB, outperforming both the global model and the nominal oracle-region pole assignment. This is the first successful local adaptation mechanism: sparse innovation proposals supply the missing spatial prior that held-out global rays alone could not provide.

The robustness sweep apps\_industrial\_breakthrough/exponential\_spline\_fri\_region\_stress\_sweep.py repeats the experiment over three projection-angle budgets and three additive noise levels. To keep the sweep diagnostic rather than combinatorial, it uses a compact region-basis candidate set containing the physical oracle family, the single-case selected family, a smoother support-$4$ family, and the best global family; the exhaustive $5^4$ region-basis search remains available as an optional mode. Across all nine stress cases, the FRI-region model improves over the best global held-out basis. The gain increases with measurement density, from a mean $+1.1667$ dB at $12$ angles to $+5.0277$ dB at $24$ angles, while the recovered region accuracy rises from $89.50\%$ to $100.00\%$. This confirms that the sparse-innovation step is not a one-off artifact: as the inverse problem receives enough projections to identify the transition set, local operator adaptation becomes reliably beneficial.

We then deliberately break the axis-aligned assumption with apps\_industrial\_breakthrough/exponential\_spline\_fri\_curve\_adaptation.py. The target field is generated by three smooth transition curves: two vertical pole-boundary curves and one top-interface curve. A row/column FRI-style derivative tracker fits low-order curve proposals from the global proxy and refines them by held-out rays. This harder test exposes the current bottleneck. The global held-out basis reaches $32.2114$ dB, while the true curved local operator assignment reaches $41.5032$ dB, proving that curved local operators have large headroom. However, the detected curved selector reaches only $31.1969$ dB despite $91.00\%$ region-label agreement and boundary RMSE $1.4606$. The failure is not the absence of local operator advantage; it is the scoring layer. Under curved imperfect labels, the held-out ray residual prefers smoother surrogate pole assignments rather than the physical local pole map. The next algorithmic step is therefore joint geometry–basis–coefficient refinement, or a region-contrastive validation score that prevents the local operator assignment from collapsing to a globally smooth surrogate.

The joint refinement runner apps\_industrial\_breakthrough/exponential\_spline\_fri\_curve\_joint\_refinement.py implements this next correction. It expands the curve-offset lattice near the best FRI proposal, fits coefficients for each candidate local operator assignment, and augments the held-out residual with an edge-consistency contrast term that rewards reconstructions whose gradient energy concentrates on the proposed sparse transition curves. This converts the curved selector from a negative result into a partial recovery: the selected joint model reaches $34.7216$ dB, a $+2.5102$ dB gain over the global held-out basis. The best candidate present in the searched family reaches $35.8310$ dB, while the true-curve oracle remains at $41.5032$ dB. Thus the scoring fix is directionally correct but incomplete. The remaining gap now separates two effects: curve localization error and the limited pole-assignment candidate family. This is the cleanest current target for further algorithmic work.

Finally, apps\_industrial\_breakthrough/exponential\_spline\_geometry\_reproduction.py isolates the more geometric point raised by the exponential-spline curve and surface literature [delgadogonzalo2012exponential]. The target is a closed harmonic curve with modes up to order three, represented from only twelve parameter samples and eight control points. The matched compact harmonic E-spline uses the pole set $\{0,\pm j2\pi/M,\pm j4\pi/M,\pm j6\pi/M\}$ and reaches dense curve RMSE $5.4457\times10^{-6}$. With the same number of control points and the same samples, a generic cubic polynomial spline reaches only $1.0256\times10^{-2}$ RMSE, and a piecewise-linear polygon reaches $3.8693\times10^{-2}$ RMSE. This is a small but important result: if the expected geometry is known to be harmonic, elliptic, spherical, cylindrical, or otherwise parametrizable by a known exponential-polynomial family, then OSNR should place that family directly in the geometric span rather than recover it indirectly through a generic volumetric grid. For NeRF-like scenes this suggests a patch-based route: segment or initialize object surfaces with geometry-reproducing parametric E-splines, then fit texture/radiance on those surfaces. For PDE domains it suggests an even cleaner route: represent both boundary geometry and the solution field in operator-matched spline spaces.

The follow-up optimizer apps\_industrial\_breakthrough/eggroll\_spline\_shape\_optimizer.py tests whether this geometric advantage can be used when the target shape is not known in closed form. Motivated by Schmitter and Unser's continuous-domain shape projectors and functional PCA construction [schmitter2018landmark], the experiment represents a closed spline curve by $16$ control points but restricts the learned geometric search to an $8$-dimensional continuous shape subspace. This is the analogue of replacing an arbitrary coordinate-field parameter vector by a learned spline-shape chart. The stochastic search component is motivated by the EGGROLL low-rank evolution-strategy result [sarkar2026eggroll]: rank-one perturbations can be evaluated as hardware-friendly low-rank updates, but the experiment separates this hardware trick from the geometric prior itself.

The result is deliberately diagnostic. Full Gaussian ES over all $32$ control coordinates reaches dense RMSE $2.9097\times10^{-2}$ after $67{,}200$ forward evaluations, while rank-one EGGROLL-style perturbations applied directly to the raw control matrix reach $3.3517\times10^{-2}$. Thus low-rank noise alone does not solve geometry discovery. When the same evaluation budget is spent inside the Schmitter-style spline subspace, the error drops to $6.9268\times10^{-3}$; rank-one EGGROLL perturbations inside that subspace reach a comparable $7.4862\times10^{-3}$. The algebraic continuous-subspace projection oracle reaches $3.5689\times10^{-4}$ in $0.293$ ms, exposing the remaining optimization gap. The practical implication is precise: the promising route is not blind evolution over arbitrary OSNR coefficients, but variable projection. Use stochastic low-rank search only for nonlinear geometry, visibility, and knot variables; solve the linear radiance or texture coefficients algebraically once a candidate geometry is proposed.

The next controlled runner, apps\_industrial\_breakthrough/eggroll\_adjoint\_variable\_projection.py, implements that variable-projection step explicitly. A one-dimensional spline boundary partitions a $40\times40$ radiance field into two continuous regions. For each candidate boundary, the code constructs the ray operator $H(\varphi)$ by projecting masked smooth atoms, eliminates the linear radiance coefficients by the adjoint normal equation $$c^\star(\varphi)=\left(H(\varphi)^\top H(\varphi)+\lambda I\right)^{-1}H(\varphi)^\top y,$$ and scores the resulting field on acquisition directions not used in the coefficient solve. This is the minimal inverse-problem analogue of a NeRF geometry/radiance separation: nonlinear geometry is searched, while linear radiance is solved in closed form.

The experiment clarifies both the opportunity and the bottleneck. The mean-geometry variable-projection baseline reaches $16.8968$ dB on held-out projections. Full Gaussian ES over raw boundary controls improves to $18.9082$ dB, rank-one EGGROLL over raw controls reaches $19.9034$ dB, and rank-one EGGROLL in the six-dimensional spline subspace reaches $20.1043$ dB. However, the true-geometry variable-projection oracle reaches $33.1874$ dB with a field RMSE of $3.3521\times10^{-2}$. Thus the adjoint variable-projection mechanism is working, but stochastic boundary discovery remains underidentified from the current projection residual alone. This is an important negative constraint for the NeRF/SIREN campaign: the next improvement must add stronger geometry evidence, such as FRI edge measurements, silhouette consistency, epipolar visibility constraints, or a learned continuous shape prior. More low-rank perturbation budget alone is unlikely to close the oracle gap.

The Haouchat-matched follow-up apps\_industrial\_breakthrough/haouchat\_matched\_variable\_projection.py performs that correction. Instead of using a primitive projection mask, it builds the inner ray operator from the same quadrature-evaluated tensor-product spline basis used in the $52$–$56$ dB matched-ray experiments. The target data are generated by a damped-harmonic exponential spline ray operator. Candidate geometries still define a boundary-dependent masked coefficient dictionary, but the forward map is now $$H_\varphi = H_{\beta_L}\,\Phi(\varphi),$$ where $H_{\beta_L}$ is the Haouchat-style ray projector for the selected tensor-product basis and $\Phi(\varphi)$ applies the boundary-dependent coefficient atoms. The stochastic outer loop is also given a weak FRI-like edge observation of the boundary, so the score combines held-in projection residual and edge consistency. This is the first experiment in this line that combines all three ingredients: matched spline rays, adjoint variable projection, and sparse boundary evidence.

The result changes the interpretation sharply. With the matched damped-harmonic operator, the true-geometry projection oracle reaches $96.9888$ dB on held-out rays, confirming that the inner ray/inverse model itself is not the limiting factor. The edge-aware rank-one spline-subspace EGGROLL search reaches $42.2339$ dB, up from the mean-geometry baseline of $33.3709$ dB, with boundary RMSE reduced from $9.5169\times10^{-2}$ to $2.7066\times10^{-2}$. The exponential-decay candidate also benefits from edge evidence, improving from $40.6125$ dB without the edge term to $41.0644$ dB with it. This supports the current thesis: OSNR does not need a dense NeRF-style MLP to represent the radiance once the operator is matched; the hard remaining problem is physically constrained geometry and visibility discovery.

The follow-up apps\_industrial\_breakthrough/haouchat\_fri\_edge\_variable\_projection.py removes the remaining artificial part of the edge-aware score. Instead of injecting a boundary hint directly from the true geometry, it forms a measurement-derived FRI proxy: training rays are backprojected through the matched adjoint, row-wise derivatives of the normalized adjoint image are localized, and the resulting peak track is smoothed into a candidate boundary. This is not a complete multidimensional FRI surface solver, but it is a measurement-only sparse-transition proposal. In the default run, the synthetic edge hint has boundary RMSE $1.3676\times10^{-2}$, while the adjoint-derived edge track has RMSE $1.9117\times10^{-2}$.

Using this measurement-derived edge evidence, the matched damped-harmonic model reaches $41.5632$ dB on held-out rays, compared with $33.3709$ dB for mean geometry and $96.9888$ dB for the true-geometry oracle. The cleaner synthetic edge hint reaches $45.9174$ dB in the same runner. The gap between $41.56$ dB and $45.92$ dB is useful: it quantifies the price of deriving geometry evidence from measurements rather than providing it externally. The experiment therefore validates the direction without hiding the remaining work. The next real-scene version should replace row-wise adjoint peaks by multi-view epipolar FRI proposals and visibility-aware surface clustering.

The next controlled refinement, apps\_industrial\_breakthrough/haouchat\_multiray\_fri\_edge\_projection.py, replaces the row-wise adjoint peak estimate by a multi-ray residual selection loop. The adjoint boundary is used only as initialization. Each boundary control point is perturbed over a small local offset lattice, and candidates are scored by the held-in variable-projection residual plus curvature and anchor penalties: $$\mathcal{S}(\varphi)= \|H_{\beta_L}\Phi(\varphi)c^\star(\varphi)-y_{\mathrm{score}}\|_2^2 +\eta\|\Delta^2\varphi\|_2^2 +\rho\|\varphi-\varphi_{\mathrm{adj}}\|_2^2.$$ This converts the crude measurement edge into a ray-consistent FRI boundary proposal without accessing the hidden target geometry.

The improvement is large. The row/adjoint edge has boundary RMSE $1.9117\times10^{-2}$; the multi-ray refined edge has RMSE $1.4848\times10^{-3}$, better than the noisy synthetic edge control ($1.3676\times10^{-2}$). With the matched damped-harmonic operator, direct variable projection on the multi-ray FRI boundary reaches $66.1488$ dB on held-out rays and field RMSE $2.5205\times10^{-3}$. The true-geometry oracle remains $96.9888$ dB, so the experiment is still controlled rather than a final SOTA benchmark, but it shows that the geometry bottleneck can be attacked algebraically by ray-consistent sparse innovation refinement. Interestingly, applying the stochastic EGGROLL search after this refined geometry is worse ($40.3169$ dB), so the current best path is not more random search but better deterministic boundary proposal.

The robustness sweep apps\_industrial\_breakthrough/haouchat\_multiray\_fri\_robustness\_sweep.py repeats the deterministic part of the experiment across two boundary variants, projection-angle counts $\{12,18,24\}$, and additive noise levels $\{0,3\times10^{-4},10^{-3}\}$. The mean-geometry baseline averages $32.5922$ dB over the $18$ cases; direct row/adjoint edge projection averages $38.0604$ dB; multi-ray FRI variable projection averages $64.0770$ dB; and the true-geometry oracle averages $97.9883$ dB. The multi-ray method remains above $61.76$ dB in every tested case. This confirms that the $66.15$ dB result is not a one-off numerical accident, but a stable consequence of selecting boundary offsets by held-in ray residuals.