Robust data-driven discovery of fractional differential equations via weak formulations and Pareto-based subset selection

# Robust data-driven discovery of fractional differential equations via weak formulations and Pareto-based subset selection
Source: [https://arxiv.org/html/2608.12879](https://arxiv.org/html/2608.12879)
Fractional partial differential equations describe nonlocal dynamics, but discovering them from noisy data is difficult because fractional differentiation amplifies high\-frequency measurement noise and the derivative orders are unknown\. We propose Weak\-Pareto, which combines an adjoint\-consistent weak formulation of fractional terms with Pareto\-based subset selection over discrete term types and continuous fractional orders\. For linear right\-hand\-side terms, the adjoint transfers fractional operators from measured fields to smooth test functions, replacing noise\-sensitive pointwise differentiation with smoothing integration; for nonlinear terms, the noise\-suppression effect is partial yet useful\. Coefficients are fitted by ridge regression within a branch\-aware differential\-evolution search over the orders\. The support size is then selected at the validation\-error–complexity elbow\. We show that the variance of fixed linear right\-hand\-side weak features vanishes under grid refinement, whereas noise amplification in pointwise fractional features increases with derivative order\. Across fractional advection–diffusion, reaction–diffusion, and Burgers benchmarks, Weak\-Pareto recovers parsimonious structures from clean and noisy measurements\. In controlled advection–diffusion and Burgers comparisons, it retains the correct support at every tested multiplicative\-noise level, whereas the unregularised strong\-form counterpart largely fails once noise is introduced; this advantage persists under additive Gaussian noise\. Ablations show that the weak library drives noise robustness and that continuous\-order Pareto search avoids the support\-selection failure of a dense fixed dictionary\. On the advection–diffusion benchmark, Weak\-Pareto yields more consistent operator recovery and substantially lower measured runtime than a neural baseline\. A two\-dimensional extension of Weak\-Pareto can recover distinct coordinate\-dependent orders, including a three\-term anisotropic model under noisy conditions\.

MSC Classification\]35R11, 65M32, 62J07, 35R30

Pongpisit ThanasutivesEmail:[pongpisit\.thanasutives@riken\.jp](mailto:)Affiliation:Center for Advanced Intelligence Project \(AIP\), RIKEN, Tokyo, JapanYoshinobu KawaharaAffiliation:Center for Advanced Intelligence Project \(AIP\), RIKEN, Tokyo, JapanAffiliation:Graduate School of Information Science and Technology, The University of Osaka, Osaka, Japan

###### keywords

fractional differential equations, data\-driven equation discovery, weak formulation, Pareto\-based subset selection, differential evolution, model selection

### Graphical Abstract

![[Uncaptioned image]](https://arxiv.org/html/2608.12879v1/GraphicalAbstract.png)

### 1Introduction

Data\-driven equation discovery aims to infer an interpretable, closed\-form model directly from observations, rather than deriving it purely from first principles\. Methods such as sparse identification of nonlinear dynamics \(SINDy\) and PDE functional identification \(PDE\-FIND\) select a small set of active terms from an overcomplete candidate library by combining regression with sparsity\-promoting model selection[1](https://arxiv.org/html/2608.12879#bib.bib1);[14](https://arxiv.org/html/2608.12879#bib.bib2);[17](https://arxiv.org/html/2608.12879#bib.bib3)\. Their central principle is parsimony: the governing equation should contain only the important terms needed to explain the observed dynamics[7](https://arxiv.org/html/2608.12879#bib.bib5)\.

Many important real\-world systems are nonlocal\. Fractional differential equations, including fractional partial differential equations \(FPDEs\), describe anomalous diffusion, viscoelasticity, transport in heterogeneous media, and several biological and financial processes through real\-valued derivative orders and nonlocal memory or spatial action[12](https://arxiv.org/html/2608.12879#bib.bib9);[10](https://arxiv.org/html/2608.12879#bib.bib11);[6](https://arxiv.org/html/2608.12879#bib.bib10)\. Discovering these equations would allow the data to reveal both the governing terms and the degree of nonlocality\.

However, fractional equation discovery introduces two difficulties beyond the integer\-order setting\. First, fractional differentiation can amplify high\-frequency measurement noise\. For the periodic spectral operators used in several benchmarks, this mechanism is explicit in the frequency domain: an operator of orderβ\\betascales a Fourier mode of wavenumberκ\\kappaby a factor whose magnitude grows as\|κ\|β\|\\kappa\|^\{\\beta\}\. Strong\-form methods therefore become increasingly fragile as the order or noise level grows\. Second, the unknown orders are continuous\. A fixed dictionary must either use a coarse order grid, which creates discretisation bias, or a dense grid, which produces many nearly collinear columns and thus destabilises support selection\.

Weak formulations address the first difficulty for integer\-order equations\. Multiplying by smooth test functions and integrating by parts transfers derivatives from the measurements to the test functions, replacing pointwise derivative estimates with integral measurements\. This principle underlies weak SINDy and related integral, Galerkin, and neural weak\-form methods, which are substantially more noise\-robust than strong\-form regression[5](https://arxiv.org/html/2608.12879#bib.bib15);[13](https://arxiv.org/html/2608.12879#bib.bib8);[8](https://arxiv.org/html/2608.12879#bib.bib6);[9](https://arxiv.org/html/2608.12879#bib.bib7);[16](https://arxiv.org/html/2608.12879#bib.bib4);[21](https://arxiv.org/html/2608.12879#bib.bib16);[3](https://arxiv.org/html/2608.12879#bib.bib17);[19](https://arxiv.org/html/2608.12879#bib.bib18)\. Complementary uncertainty\-aware model\-selection methods, including the uncertainty\-penalised Bayesian information criterion \(UBIC\), improve robustness by penalising candidate terms whose inferred coefficients exhibit high uncertainty[22](https://arxiv.org/html/2608.12879#bib.bib21)\.

Existing fractional\-discovery methods, however, still evaluate fractional derivatives pointwise, either directly on a mesh or after reconstructing the field with a neural network[4](https://arxiv.org/html/2608.12879#bib.bib12);[11](https://arxiv.org/html/2608.12879#bib.bib13);[23](https://arxiv.org/html/2608.12879#bib.bib19)\. Extending the weak formulation to fractional operators is not trivial, because Caputo, Riemann–Liouville, Grünwald–Letnikov, Riesz, and periodic spectral derivatives have different adjoints and boundary terms\. A weak feature is valid only when it uses the adjoint of the operator being identified\. To address the aforementioned difficulties, we propose*Weak\-Pareto*, a data\-driven discovery framework that couples an adjoint\-consistent weak formulation of the fractional candidate library with continuous\-order, Pareto\-based subset selection\. Table[1](https://arxiv.org/html/2608.12879#S1.T1)characterises the differences between Weak\-Pareto and existing equation\-discovery methods\. Note that notable symbolic\-regression discovery frameworks[18](https://arxiv.org/html/2608.12879#bib.bib22);[2](https://arxiv.org/html/2608.12879#bib.bib23)search over algebraic expression trees built from a fixed operator set, and usually do not support a continuously varying differential order\.

Table 1:Comparison with prior equation\-discovery methods\. “Cont\. order” means real\-valued order optimisation rather than selection from a fixed order dictionary\. “Search strategy” describes how active terms or a prescribed model are identified; “Frac\. adjoint” means that weak features use the adjoint and boundary terms of the represented fractional operator\.To our knowledge, Weak\-Pareto is the first equation\-discovery framework to achieve all five properties in Table[1](https://arxiv.org/html/2608.12879#S1.T1)\. Our core contributions are threefold\.

1. 1\.An adjoint\-consistent weak library for fractional operators\(Section[3](https://arxiv.org/html/2608.12879#S3)\)\. Each candidate uses the adjoint and boundary terms of its declared operator\. For linear right\-hand\-side terms, the fractional derivative acts only on smooth test functions; those library columns therefore do not differentiate the measured field\. Proposition[1](https://arxiv.org/html/2608.12879#Thmtheorem1)shows that the variance of a fixed linear right\-hand\-side weak feature vanishes under grid refinement, whereas the variance of a pointwise positive\-order feature diverges at a rate that increases with derivative order\. For nonlinear terms, we show explicitly that the weak row is an averaged projection of the corresponding strong feature and quantify the resulting noise\-induced bias\. For a fixed support, Proposition[3](https://arxiv.org/html/2608.12879#Thmtheorem3)further characterises when a spatial\-order perturbation can be absorbed, to first order, by refitting the active coefficients\.
2. 2\.A branch\-aware, continuous\-order Pareto search\(Section[4](https://arxiv.org/html/2608.12879#S4)\)\. Candidate terms are encoded as integer powers paired with continuous fractional orders\. Coefficients are fitted analytically inside a differential\-evolution search over the orders, while subunit, exact\-integer, and superunit Caputo modes are treated as distinct branches\. The support size is then selected at the validation\-error–complexity elbow\. This avoids both the discretisation error imposed by coarse fixed dictionaries and the severe collinearity of dense ones\. After selecting the support and temporal branch, we evaluate the operators at locally refined orders and refit the coefficients on the full data\. This removes interpolation error from the reported coefficients, although the selected model can still depend on the original order grid and search trajectory\.
3. 3\.Controlled analyses and baseline comparisons\(Section[5](https://arxiv.org/html/2608.12879#S5)\)\. Library\-only ablations isolate the advantage of weak over pointwise measurements; fixed\-dictionary ablations isolate the advantage of continuous\-order search; and a controlled comparison assesses an adapted neural fractional\-discovery framework on the advection–diffusion benchmark\. A semi\-analytic fixed\-support diagnostic tests the previously unreported superunit temporal branch\. The experiments also identify the method’s present limits on Riesz reaction–diffusion problems, demonstrate weak\-form integral fitting on irregular frozen\-soil creep data, and use an anisotropic two\-dimensional example to showcase direction\-labelled continuous\-order encoding\.

On the fractional advection–diffusion \(FADE\) and fractional Burgers benchmarks, Weak\-Pareto recovers the correct support in all five seeds at every tested multiplicative\-noise level up to20%20\\%\. Under matched selection, the strong\-form framework achieves no complete operator recovery in any noisy run\. The two\-dimensional example recovers the support, coordinate directions, and derivative orders in all 25 runs for each of its two benchmarks\.

Fig\.[1](https://arxiv.org/html/2608.12879#S1.F1)connects the main stages of Weak\-Pareto to concrete results from the supplied experiments\. Panel \(a\) visualises the FADE benchmark under10%10\\%additive Gaussian noise\. Panel \(b\) combines the adjoint transfer used for linear weak features with the branch\-aware Pareto search over support size, discrete term types, and continuous fractional orders\. Panel \(c\) displays the closed\-form encoding, the correct\-support recovery pattern in the additive\-Gaussian comparisons, the representative FADE discovery reported in Appendix[15](https://arxiv.org/html/2608.12879#S15), and an example noisy estimate for the three\-term anisotropic two\-dimensional benchmark \([28](https://arxiv.org/html/2608.12879#S5.E28)\)\. The recovery counts refer specifically to correct\-support recovery; complete operator recovery additionally requires the fractional orders to satisfy the declared tolerances\.

![Refer to caption](https://arxiv.org/html/2608.12879v1/Fig1.png)Figure 1:Overview of Weak\-Pareto\. \(a\) Visualisation of the FADE benchmark with10%10\\%additive Gaussian noise, illustrating the high\-frequency perturbations amplified by fractional differentiation\. \(b\) For a linear operator𝒳β\\mathcal\{X\}\_\{\\beta\}, the adjoint identity transfers differentiation from the measured field to a smooth test function\. The lower schematic shows the branch\-aware Pareto search over support size, discrete term types, and continuous fractional orders, with the selected validation\-error–complexity elbow highlighted\. Nonlinear weak features remain averaged measurements; therefore, their noise suppression is partial\. \(c\) The candidate encoding and its discoveries remain explicit closed\-form equations\. For both the FADE and fractional Burgers experiments at10%10\\%additive Gaussian noise, Weak\-Pareto recovers the correct support in5/55/5runs, whereas the strong\-form counterpart fails in all five runs \(Appendix[13](https://arxiv.org/html/2608.12879#S13)\)\. By augmenting the encoding tuple with a direction label, we extend Weak\-Pareto to facilitate data\-driven FPDE discovery in two\-dimensional settings \([28](https://arxiv.org/html/2608.12879#S5.E28)\)The remainder of the paper defines the model class \(Section[2](https://arxiv.org/html/2608.12879#S2)\), develops the weak library and its noise analysis \(Section[3](https://arxiv.org/html/2608.12879#S3)\), presents the Pareto search \(Section[4](https://arxiv.org/html/2608.12879#S4)\), reports the experiments \(Section[5](https://arxiv.org/html/2608.12879#S5)\), and concludes in Section[6](https://arxiv.org/html/2608.12879#S6)\. Proofs, numerical verification, and detailed settings are provided in the appendices\.

### 2Problem formulation

#### 2\.1Fractional derivative operators

We first consider a scalar fieldu⁡\(t,x\)u\(t,x\)onQ=\(0,T\)×ΩQ=\(0,T\)\\times\\Omega, withΩ⊂ℝ\\Omega\\subset\\mathbb\{R\}, observed on a space–time grid\. This notation covers the core method and analysis; Section[5\.8](https://arxiv.org/html/2608.12879#S5.SS8)gives an example extension to two spatial coordinates\. The temporal operator is denoted by𝒯mα,α\\mathcal\{T\}\_\{m\_\{\\alpha\},\\alpha\}, wheremαm\_\{\\alpha\}identifies the Caputo branch andα\\alphaits order\. We distinguish the subunit branchmα=subm\_\{\\alpha\}=\\mathrm\{sub\}for0<α<10<\\alpha<1, the exact integer modemα=intm\_\{\\alpha\}=\\mathrm\{int\}atα=1\\alpha=1, and, when the declared range extends above one, the superunit branchmα=supm\_\{\\alpha\}=\\mathrm\{sup\}for1<α<21<\\alpha<2\. Spatial terms use operators𝒳β\\mathcal\{X\}\_\{\\beta\}of orderβ\\beta\. We adopt the following standard definitions[12](https://arxiv.org/html/2608.12879#bib.bib9);[6](https://arxiv.org/html/2608.12879#bib.bib10)\. Forμ\>0\\mu\>0and nonintegerγ\>0\\gamma\>0, letn=⌈γ⌉n=\\lceil\\gamma\\rceil\. The left and right Riemann–Liouville \(RL\) integrals and derivatives on\[a,b\]\[a,b\]are

\(Iμza​f\)​\(z\)\\displaystyle\(\{\}\_\{a\}I\_\{z\}^\{\\mu\}f\)\(z\)=1Γ⁡\(μ\)​∫az\(z−s\)μ−1​f​\(s\)​ds,\\displaystyle=\\frac\{1\}\{\\Gamma\(\\mu\)\}\\int\_\{a\}^\{z\}\(z\-s\)^\{\\mu\-1\}f\(s\)\\,\\mathrm\{d\}s,\(Iμbz​f\)​\(z\)\\displaystyle\(\{\}\_\{z\}I\_\{b\}^\{\\mu\}f\)\(z\)=1Γ⁡\(μ\)​∫zb\(s−z\)μ−1​f​\(s\)​ds,\\displaystyle=\\frac\{1\}\{\\Gamma\(\\mu\)\}\\int\_\{z\}^\{b\}\(s\-z\)^\{\\mu\-1\}f\(s\)\\,\\mathrm\{d\}s,\(1\)\(Dγza​f\)​\(z\)\\displaystyle\(\{\}\_\{a\}D\_\{z\}^\{\\gamma\}f\)\(z\)=dnd​zn​\(In−γza​f\)​\(z\),\\displaystyle=\\frac\{\\,\\mathrm\{d\}^\{n\}\}\{\\,\\mathrm\{d\}z^\{n\}\}\\,\(\{\}\_\{a\}I\_\{z\}^\{n\-\\gamma\}f\)\(z\),\(Dγbz​f\)​\(z\)\\displaystyle\(\{\}\_\{z\}D\_\{b\}^\{\\gamma\}f\)\(z\)=\(−1\)n​dnd​zn​\(In−γbz​f\)​\(z\)\.\\displaystyle=\(\-1\)^\{n\}\\frac\{\\,\\mathrm\{d\}^\{n\}\}\{\\,\\mathrm\{d\}z^\{n\}\}\\,\(\{\}\_\{z\}I\_\{b\}^\{n\-\\gamma\}f\)\(z\)\.HereΓ⁡\(μ\)=∫0∞yμ−1​e−y​𝑑y\\Gamma\(\\mu\)=\\int\_\{0\}^\{\\infty\}y^\{\\mu\-1\}e^\{\-y\}\\,\\mathrm\{d\}yis Euler’s gamma function; integer\-order derivatives have their usual meaning\. The Caputo derivative subtracts the initial Taylor polynomial,DzαaC​f=Dαza​\[f−Pn−1,a​f\]\{\}\_\{a\}^\{C\}\\\!D\_\{z\}^\{\\alpha\}f=\{\}\_\{a\}D\_\{z\}^\{\\alpha\}\[f\-P\_\{n\-1,a\}f\]withPn−1,a​f​\(z\)=∑q=0n−1f\(q\)​\(a\)q\!​\(z−a\)qP\_\{n\-1,a\}f\(z\)=\\sum\_\{q=0\}^\{n\-1\}\\frac\{f^\{\(q\)\}\(a\)\}\{q\!\}\(z\-a\)^\{q\}andn=⌈α⌉n=\\lceil\\alpha\\rceil\. Below,DtαD\_\{t\}^\{\\alpha\}abbreviatesDtα0C\{\}\_\{0\}^\{C\}\\\!D\_\{t\}^\{\\alpha\}for nonintegerα\\alpha, while∂t\\partial\_\{t\}denotes the exact integer operator\. The right RL derivatives in Eq\. \([1](https://arxiv.org/html/2608.12879#S2.E1)\) arise as the adjoints of the corresponding left\-sided operators \(Section[3\.2](https://arxiv.org/html/2608.12879#S3.SS2)\)\. On periodic domains we also use the Riesz operator and the directional spectral derivative, defined through their Fourier multipliers

ℱ⁡\[ℛβ​f\]​\(κ\)=−\|κ\|β​ℱ​\[f\]​\(κ\),ℱ⁡\[Dxβ​f\]​\(κ\)=\(i​κ\)β​ℱ​\[f\]​\(κ\)\.\\mathcal\{F\}\[\\mathcal\{R\}\_\{\\beta\}f\]\(\\kappa\)=\-\|\\kappa\|^\{\\beta\}\\,\\mathcal\{F\}\[f\]\(\\kappa\),\\hskip 20\.00003pt\\mathcal\{F\}\[D\_\{x\}^\{\\beta\}f\]\(\\kappa\)=\(\\mathrm\{i\}\\kappa\)^\{\\beta\}\\,\\mathcal\{F\}\[f\]\(\\kappa\)\.\(2\)Hereℱ\\mathcal\{F\}denotes the Fourier transform andκ\\kappais the spatial wavenumber\. The Riesz operatorℛβ=−\(−Δ\)β/2\\mathcal\{R\}\_\{\\beta\}=\-\(\-\\Delta\)^\{\\beta/2\}has the non\-positive multiplier−\|κ\|β\-\|\\kappa\|^\{\\beta\}; thusβ=2\\beta=2recovers∂x2\\partial\_\{x\}^\{2\}, and a positive coefficient produces diffusion\. The directional multiplier\(i​κ\)β\(\\mathrm\{i\}\\kappa\)^\{\\beta\}uses the principal branch,

\(i​κ\)β:=\|κ\|β​exp⁡\(i​π2​β​sgn⁡\(κ\)\),κ≠0,\(\\mathrm\{i\}\\kappa\)^\{\\beta\}:=\|\\kappa\|^\{\\beta\}\\exp\\\!\\Bigl\(\\mathrm\{i\}\\tfrac\{\\pi\}\{2\}\\beta\\,\\operatorname\{sgn\}\(\\kappa\)\\Bigr\),\\hskip 20\.00003pt\\kappa\\neq 0,\(3\)wheresgn\\operatorname\{sgn\}is the sign function and the zero mode is annihilated\. The multiplier is conjugate\-symmetric,\(i​κ\)β¯=\(i⁡\(−κ\)\)β\\overline\{\(\\mathrm\{i\}\\kappa\)^\{\\beta\}\}=\(\\mathrm\{i\}\(\-\\kappa\)\)^\{\\beta\}; consequently, it maps real fields to real fields and models advective fractional transport;β=1\\beta=1recovers∂x\\partial\_\{x\}\.

The searched model class is the parsimonious fractional equation

𝒯mα,α​u=∑j=1cξj​upj​𝒳βj​u\.\\mathcal\{T\}\_\{m\_\{\\alpha\},\\alpha\}u=\\sum\_\{j=1\}^\{c\}\\xi\_\{j\}\\,u^\{p\_\{j\}\}\\,\\mathcal\{X\}\_\{\\beta\_\{j\}\}u\.\(4\)The model hasccactive terms, with integer powerspjp\_\{j\}, spatial ordersβj\\beta\_\{j\}, and coefficientsξj\\xi\_\{j\}\. The class includes fractional advection–diffusion \(pj=0p\_\{j\}=0\) and nonlinear transport such as the Burgers termu​∂xuu\\,\\partial\_\{x\}u\(p=1,β=1p=1,\\beta=1\)\. We define𝒳0\\mathcal\{X\}\_\{0\}as the identity; hence\(p,0\)\(p,0\)represents the reaction termup\+1u^\{p\+1\}\. This is a modelling convention, not theβ→0\\beta\\to 0limit of the Riesz multiplier, which is−\(Id−Π0\)\-\(\\mathrm\{Id\}\-\\Pi\_\{0\}\)because the zero Fourier mode is annihilated; it reduces to−Id\-\\mathrm\{Id\}only on mean\-zero fields, and any sign is absorbed intoξj\\xi\_\{j\}\.

#### 2\.2FPDE encoding

Each candidate is represented by the tuple

ℳ=\(mα,α,𝒑,𝜷,𝝃\),𝒑=\(p1,…,pc\),𝜷=\(β1,…,βc\),𝝃=\(ξ1,…,ξc\)\.\\mathcal\{M\}=\(m\_\{\\alpha\},\\alpha,\\bm\{p\},\\bm\{\\beta\},\\bm\{\\xi\}\),\\hskip 20\.00003pt\\bm\{p\}=\(p\_\{1\},\\dots,p\_\{c\}\),\\hskip 10\.00002pt\\bm\{\\beta\}=\(\\beta\_\{1\},\\dots,\\beta\_\{c\}\),\\hskip 10\.00002pt\\bm\{\\xi\}=\(\\xi\_\{1\},\\dots,\\xi\_\{c\}\)\.\(5\)Heremαm\_\{\\alpha\}records the temporal branch\. This label is crucial at an integer: an estimate such asα=0\.999\\alpha=0\.999remains a subunit Caputo model and is not identified with∂t\\partial\_\{t\}\. The subunit, exact\-integer, and superunit branch minima are therefore optimised separately and compared using the same validation objective\.

Each right\-hand\-side term is represented by a discrete powerpjp\_\{j\}and a continuous orderβj\\beta\_\{j\}\. Unlike fixed\-library regression, this encoding does not enumerate all possible orders in advance\. A finite fixed grid either omits the true order or becomes increasingly collinear as it is refined \(Section[5\.10](https://arxiv.org/html/2608.12879#S5.SS10)\); Weak\-Pareto instead assembles only the columns proposed by the search\. For example, the FADE equationDt0\.8u=−∂xu\+0\.5Dx1\.7uD\_\{t\}^\{0\.8\}u=\-\\partial\_\{x\}u\+0\.5D\_\{x\}^\{1\.7\}uis encoded bymα=subm\_\{\\alpha\}=\\mathrm\{sub\},α=0\.8\\alpha=0\.8,𝒑=\(0,0\)\\bm\{p\}=\(0,0\),𝜷=\(1,1\.7\)\\bm\{\\beta\}=\(1,1\.7\), and𝝃=\(−1,0\.5\)\\bm\{\\xi\}=\(\-1,0\.5\)\. At each support sizecc, the method searches for the bestcc\-term equation\. Within a fractional temporal branch it optimises\(α,β1,…,βc\)\(\\alpha,\\beta\_\{1\},\\dots,\\beta\_\{c\}\); in the exact integer mode,α=1\\alpha=1is fixed and only𝜷\\bm\{\\beta\}is optimised\. For fixed modes, powers, and orders, the coefficients follow from linear regression \(Section[4](https://arxiv.org/html/2608.12879#S4)\)\. In multiple spatial coordinates, our encoding can be augmented by a direction labeldjd\_\{j\}for each term\. The example in Section[5\.8](https://arxiv.org/html/2608.12879#S5.SS8)usesdj∈\{x,y\}d\_\{j\}\\in\\\{x,y\\\}and searches the direction pattern together with the continuous orders; its reported implementation is restricted to linear terms\.

### 3Weak formulation of fractional candidate terms

This section develops the first core contribution: weak features that are consistent with the adjoint and boundary structure of each fractional operator, together with a feature\-level explanation of their noise\-robustness advantage over pointwise differentiation\.

#### 3\.1Test functions and weak features

Let\{ϕk\}k=1K\\\{\\phi\_\{k\}\\\}\_\{k=1\}^\{K\}be smooth, separable test functions

ϕk​\(t,x\)=ϑℓt​\(t\)​ψℓx​\(x\),ϕk∈Cs​\(Q¯\)\.\\phi\_\{k\}\(t,x\)=\\vartheta\_\{\\ell\_\{t\}\}\(t\)\\,\\psi\_\{\\ell\_\{x\}\}\(x\),\\hskip 20\.00003pt\\phi\_\{k\}\\in C^\{s\}\(\\overline\{Q\}\)\.\(6\)The indexkkenumerates temporal–spatial test\-function pairs\(ℓt,ℓx\)\(\\ell\_\{t\},\\ell\_\{x\}\)\. For a localised test family, we use the term*test window*for a test function together with the region on which it has appreciable weight\. Heres∈ℕ0s\\in\\mathbb\{N\}\_\{0\},Q¯=\[0,T\]×Ω¯\\overline\{Q\}=\[0,T\]\\times\\overline\{\\Omega\}, andCs​\(Q¯\)C^\{s\}\(\\overline\{Q\}\)denotes functions whose partial derivatives of total order at mostssextend continuously toQ¯\\overline\{Q\}; the smoothness order is chosen for the highest derivative in the library\. The implementation provides compactly supported bumps, localised Gaussian test windows, and global Fourier modes\. Reported experiments use Gaussian test windows and, for periodic high\-order Riesz cases, Fourier spatial modes; bumps are not used in the reported results\. In this subsection, “compact support” refers only to localisation of a test function\. Elsewhere, the support of a candidate model means its active\-term set in the subset search\. Compact support of a test function or periodicity removes the corresponding continuum boundary terms\. Gaussian tests use the exact discrete adjoint of Section[3\.3](https://arxiv.org/html/2608.12879#S3.SS3), preserving the discrete identity without a zero\-trace assumption\. Each test function produces one scalar equation,

⟨𝒯mα,α​u,ϕk⟩=∑j=1cξj​⟨upj​𝒳βj​u,ϕk⟩,⟨f,g⟩=∫Qf​g​𝑑x​𝑑t\.\\langle\\mathcal\{T\}\_\{m\_\{\\alpha\},\\alpha\}u,\\,\\phi\_\{k\}\\rangle=\\sum\_\{j=1\}^\{c\}\\xi\_\{j\}\\,\\langle u^\{p\_\{j\}\}\\mathcal\{X\}\_\{\\beta\_\{j\}\}u,\\,\\phi\_\{k\}\\rangle,\\hskip 20\.00003pt\\langle f,\\,g\\rangle=\\int\_\{Q\}f\\,g\\,\\mathrm\{d\}x\\,\\mathrm\{d\}t\.\(7\)The construction rests on the adjoint identity

⟨ℒ​f,ϕk⟩=⟨f,ℒ∗​ϕk⟩\+Bℒ​\(f,ϕk\),\\langle\\mathcal\{L\}f,\\,\\phi\_\{k\}\\rangle=\\langle f,\\,\\mathcal\{L\}^\{\\ast\}\\phi\_\{k\}\\rangle\+B\_\{\\mathcal\{L\}\}\(f,\\phi\_\{k\}\),\(8\)whereℒ∗\\mathcal\{L\}^\{\\ast\}is the adjoint andBℒB\_\{\\mathcal\{L\}\}contains boundary and initial contributions\. The sameϕk\\phi\_\{k\}is used on both sides; the operatorℒ\\mathcal\{L\}determines its adjoint and boundary contribution\. Compact support of a test function or periodicity removes continuum spatial boundary terms\. For Gaussian tests, the exact discrete adjoint replaces a zero\-trace assumption; the Caputo target retains its initial\-condition correction\. Linear spatial operators transfer completely from the data to the test function\. Forup​𝒳β​uu^\{p\}\\mathcal\{X\}\_\{\\beta\}u, the discrete adjoint gives an integrated projection of the strong nonlinear feature, averaging it without removing differentiation from the noisy data path\.

#### 3\.2Adjoint identities

##### Temporal \(Caputo\) target\.

Letn=⌈α⌉n=\\lceil\\alpha\\rceilandPn−1,0​u​\(t,x\)=∑q=0n−1∂tqu⁡\(0,x\)​tq/q\!P\_\{n\-1,0\}u\(t,x\)=\\sum\_\{q=0\}^\{n\-1\}\\partial\_\{t\}^\{q\}u\(0,x\)t^\{q\}/q\!\. For0<α<20<\\alpha<2,α≠1\\alpha\\neq 1, and test functions satisfying the terminal conditions required to remove the right\-endpoint traces, fractional integration by parts gives

⟨Dtα0C​u,ϕk⟩=⟨u−Pn−1,0​u,\(DαTt\)​ϕk⟩\.\\langle\{\}\_\{0\}^\{C\}\\\!D\_\{t\}^\{\\alpha\}u,\\,\\phi\_\{k\}\\rangle=\\langle\\,u\-P\_\{n\-1,0\}u\\,,\\,\(\{\}\_\{t\}D\_\{T\}^\{\\alpha\}\)\\phi\_\{k\}\\rangle\.\(9\)Thus the subunit branch subtractsu⁡\(0,⋅\)u\(0,\\cdot\), whereas the superunit branch also subtractst​∂tu⁡\(0,⋅\)t\\,\\partial\_\{t\}u\(0,\\cdot\)\. Atα=1\\alpha=1, the trace\-free continuum identity is∫Qutϕk=−∫Qu∂tϕk\\int\_\{Q\}u\_\{t\}\\phi\_\{k\}=\-\\int\_\{Q\}u\\,\\partial\_\{t\}\\phi\_\{k\}; reported Gaussian tests instead use the discrete transpose identity of Section[3\.3](https://arxiv.org/html/2608.12879#S3.SS3)\. The domain integral leavesuuundifferentiated, although the superunit branch still depends on∂tu⁡\(0,⋅\)\\partial\_\{t\}u\(0,\\cdot\)\. In the reported L1 experiments,Dhα=𝖫𝟣α−1​D1D\_\{h\}^\{\\alpha\}=\\mathsf\{L1\}\_\{\\alpha\-1\}D\_\{1\}and\(Dhα\)⊤=D1⊤​𝖫𝟣α−1⊤\(D\_\{h\}^\{\\alpha\}\)^\{\\top\}=D\_\{1\}^\{\\top\}\\mathsf\{L1\}\_\{\\alpha\-1\}^\{\\top\}\. Its first weights are endpoint\-concentrated, implicitly treating the initial rate one\-sidedly\. Hence the target is weak in the interior but not derivative\-free at the initial boundary\. Remark[2](https://arxiv.org/html/2608.12879#Thmremark2)gives the scaling and noise analysis\. The optional Volterra target estimates the initial Taylor polynomial explicitly but is not used here\.

##### Spatial terms\.

For a linear spatial term the weak feature is

θk,β=⟨𝒳β​u,ϕk⟩=⟨u,𝒳β∗​ϕk⟩,\\theta\_\{k,\\beta\}=\\langle\\mathcal\{X\}\_\{\\beta\}u,\\,\\phi\_\{k\}\\rangle=\\langle u,\\,\\mathcal\{X\}\_\{\\beta\}^\{\\ast\}\\phi\_\{k\}\\rangle,\(10\)where the adjoint depends on the operator definition\. For a left Riemann–Liouville derivative inxxon\[a,b\]\[a,b\], suppressing the other variables, the continuum identity is

∫ab\(Dβxa​u\)​\(x\)​ϕ​\(x\)​𝑑x=∫abu⁡\(x\)​\(Dβbx​ϕ\)​\(x\)​𝑑x,\\int\_\{a\}^\{b\}\(\{\}\_\{a\}D\_\{x\}^\{\\beta\}u\)\(x\)\\,\\phi\(x\)\\,\\mathrm\{d\}x=\\int\_\{a\}^\{b\}u\(x\)\\,\(\{\}\_\{x\}D\_\{b\}^\{\\beta\}\\phi\)\(x\)\\,\\mathrm\{d\}x,under the endpoint conditions stated in Appendix[7](https://arxiv.org/html/2608.12879#S7); hence\(Dβxa\)∗=Dβbx\(\{\}\_\{a\}D\_\{x\}^\{\\beta\}\)^\{\\ast\}=\{\}\_\{x\}D\_\{b\}^\{\\beta\}on that domain\. For the periodic Riesz operator,ℛβ∗=ℛβ\\mathcal\{R\}\_\{\\beta\}^\{\\ast\}=\\mathcal\{R\}\_\{\\beta\}because its multiplier−\|κ\|β\-\|\\kappa\|^\{\\beta\}is real and even\. For the directional periodic operator,

ℱ⁡\[𝒳β∗​ϕk\]​\(κ\)=\(i​κ\)β¯​ℱ​\[ϕk\]​\(κ\)=\(i⁡\(−κ\)\)β​ℱ​\[ϕk\]​\(κ\)\.\\mathcal\{F\}\[\\mathcal\{X\}\_\{\\beta\}^\{\\ast\}\\phi\_\{k\}\]\(\\kappa\)=\\overline\{\(\\mathrm\{i\}\\kappa\)^\{\\beta\}\}\\,\\mathcal\{F\}\[\\phi\_\{k\}\]\(\\kappa\)=\(\\mathrm\{i\}\(\-\\kappa\)\)^\{\\beta\}\\,\\mathcal\{F\}\[\\phi\_\{k\}\]\(\\kappa\)\.Using an adjoint from a different operator family changes the model being identified and can bias the recovered order\. Each benchmark therefore uses the adjoint of its declared candidate operator\.

##### Nonlinear terms\.

For the nonlinear termup​𝒳β​uu^\{p\}\\mathcal\{X\}\_\{\\beta\}uin \([4](https://arxiv.org/html/2608.12879#S2.E4)\),

⟨up​𝒳β​u,ϕk⟩=⟨𝒳β​u,up​ϕk⟩=⟨u,𝒳β∗​\(up​ϕk\)⟩\.\\langle u^\{p\}\\mathcal\{X\}\_\{\\beta\}u,\\,\\phi\_\{k\}\\rangle=\\langle\\mathcal\{X\}\_\{\\beta\}u,\\,u^\{p\}\\phi\_\{k\}\\rangle=\\langle u,\\,\\mathcal\{X\}\_\{\\beta\}^\{\\ast\}\(u^\{p\}\\phi\_\{k\}\)\\rangle\.\(11\)Eq\. \([11](https://arxiv.org/html/2608.12879#S3.E11)\) is the consistent weak form of the nonlinear candidate\. At the discrete level, letAArepresent𝒳β\\mathcal\{X\}\_\{\\beta\}\. Eq\. \([12](https://arxiv.org/html/2608.12879#S3.E12)\) gives the exact chain

⟨u,A∗,h​\(up​ϕk\)⟩h=⟨A​u,up​ϕk⟩h=⟨up​A​u,ϕk⟩h\.\\langle u,\\,A^\{\\ast,h\}\(u^\{p\}\\phi\_\{k\}\)\\rangle\_\{h\}=\\langle Au,\\,u^\{p\}\\phi\_\{k\}\\rangle\_\{h\}=\\langle u^\{p\}Au,\\,\\phi\_\{k\}\\rangle\_\{h\}\.Thus the discrete weak nonlinear feature on the left is exactly the test\-function projection of the corresponding pointwise strong featureup​A​uu^\{p\}Au\. For nonlinear terms, differentiation still acts on the measured field, while integration provides averaging\. The resulting feature may therefore be biased, and the strongest noise guarantee applies to linear terms\.

#### 3\.3Numerical computation of adjoint identities

The numerical library uses the discrete adjoint of the operator that defines each candidate\. LetA∈ℝn×nA\\in\\mathbb\{R\}^\{n\\times n\}be the real\-grid matrix representing a one\-dimensional discrete operator, letW=diag⁡\(w0,…,wn−1\)∈ℝn×nW=\\operatorname\{diag\}\(w\_\{0\},\\ldots,w\_\{n\-1\}\)\\in\\mathbb\{R\}^\{n\\times n\}be the positive quadrature\-weight matrix, and letf,ϕ∈ℝnf,\\phi\\in\\mathbb\{R\}^\{n\}, with⟨f,g⟩h=f⊤​W​g\\langle f,\\,g\\rangle\_\{h\}=f^\{\\\!\\top\}Wg\. The discrete adjoint is

A∗,h=W−1​A⊤​W,⟨A​f,ϕ⟩h=⟨f,A∗,h​ϕ⟩h\.A^\{\\ast,h\}=W^\{\-1\}A^\{\\\!\\top\}W,\\hskip 20\.00003pt\\langle Af,\\,\\phi\\rangle\_\{h\}=\\langle f,\\,A^\{\\ast,h\}\\phi\\rangle\_\{h\}\.\(12\)WhenA=Aβ,hA=A\_\{\\beta,h\}discretises𝒳β\\mathcal\{X\}\_\{\\beta\}, we write𝒳β∗,h​ϕ:=Aβ,h∗,h​ϕ\\mathcal\{X\}\_\{\\beta\}^\{\\ast,h\}\\phi:=A\_\{\\beta,h\}^\{\\ast,h\}\\phi\. For the uniform quadrature used here,WWis a scalar multiple of the identity andA∗,h=A⊤A^\{\\ast,h\}=A^\{\\top\}\. A left Grünwald–Letnikov derivative therefore contributes\(GLγ\)⊤​ϕ\(G\_\{L\}^\{\\gamma\}\)^\{\\top\}\\phi; periodic operators use the conjugate Fourier multiplier; and the Caputo target uses the transpose of its L1 matrix\. Because the Caputo family changes definition at integer orders, the subunit branch, exact integer operator, and superunit branch are precomputed and searched separately\. Interpolation and local polishing remain within the selected branch\. The exactα=1\\alpha=1candidate is compared directly with the fractional\-branch minima, rather than approximated by a nearby fractional order\. The superunit implementation is included and verified numerically\. Appendix[8](https://arxiv.org/html/2608.12879#S8)verifies the adjoints, branch separation, and discretisation accuracy\.

#### 3\.4The noise\-robust weak library

Letu⋆u^\{\\star\}denote the noise\-free field and let the measured field beu=u⋆\+ηu=u^\{\\star\}\+\\eta, whereη\\etais zero\-mean measurement noise\. For a linear candidate \(p=0p=0\), define the noise\-free weak feature byθk,β⋆:=⟨u⋆,𝒳β∗​ϕk⟩\\theta^\{\\star\}\_\{k,\\beta\}:=\\langle u^\{\\star\},\\,\\mathcal\{X\}\_\{\\beta\}^\{\\ast\}\\phi\_\{k\}\\rangleand the measured feature byθk,β:=⟨u,𝒳β∗​ϕk⟩\\theta\_\{k,\\beta\}:=\\langle u,\\,\\mathcal\{X\}\_\{\\beta\}^\{\\ast\}\\phi\_\{k\}\\rangle\. Their difference is

θk,β−θk,β⋆=⟨η,𝒳β∗​ϕk⟩,\\theta\_\{k,\\beta\}\-\\theta^\{\\star\}\_\{k,\\beta\}=\\langle\\eta,\\,\\mathcal\{X\}\_\{\\beta\}^\{\\ast\}\\phi\_\{k\}\\rangle,\(13\)an integral projection of the noise, whereas the strong form applies the fractional operator to the noise pointwise\. Here*pointwise*means that𝒳β​u\\mathcal\{X\}\_\{\\beta\}uis evaluated at individual grid nodes before regression, rather than first being integrated or projected against a test function\. For the periodic spatial operators used below, write the Fourier action as

ℱx​\[𝒳β​f\]​\(κ\)=sβ​\(κ\)​ℱx​\[f\]​\(κ\),sβ​\(κ\)=\{−\|κ\|β,𝒳β=ℛβ,\(i​κ\)β,𝒳β=Dxβ,\\mathcal\{F\}\_\{x\}\[\\mathcal\{X\}\_\{\\beta\}f\]\(\\kappa\)=s\_\{\\beta\}\(\\kappa\)\\,\\mathcal\{F\}\_\{x\}\[f\]\(\\kappa\),\\hskip 20\.00003pts\_\{\\beta\}\(\\kappa\)=\\begin\{cases\}\-\|\\kappa\|^\{\\beta\},&\\mathcal\{X\}\_\{\\beta\}=\\mathcal\{R\}\_\{\\beta\},\\\\ \(\\mathrm\{i\}\\kappa\)^\{\\beta\},&\\mathcal\{X\}\_\{\\beta\}=D\_\{x\}^\{\\beta\},\\end\{cases\}hence both positive\-order multipliers satisfy\|sβ​\(κ\)\|=\|κ\|β\|s\_\{\\beta\}\(\\kappa\)\|=\|\\kappa\|^\{\\beta\}\. The weak–strong contrast can therefore be stated precisely for independent grid noise\.

##### Theoretical properties\.

The analysis below separates three properties\. First, adjoint consistency ensures that every weak column represents the declared fractional operator and its boundary convention\. Second, for linear right\-hand\-side columns, integration averages independent measurement noise: Proposition[1](https://arxiv.org/html/2608.12879#Thmtheorem1)shows variance decay under grid refinement, while the corresponding unregularised strong feature becomes increasingly noisy\. These results characterise individual features; estimator\-level behaviour is assessed empirically\. Corollary[2](https://arxiv.org/html/2608.12879#Thmtheorem2)extends the variance result to the principal multiplicative\-noise model, while Remarks[1](https://arxiv.org/html/2608.12879#Thmremark1)and[2](https://arxiv.org/html/2608.12879#Thmremark2)describe the remaining nonlinear bias and Caputo endpoint sensitivity\.

###### Proposition 1\.

Letηi​j\\eta\_\{ij\},i=0,…,nt−1i=0,\\dots,n\_\{t\}\-1andj=0,…,nx−1j=0,\\dots,n\_\{x\}\-1, be independent, zero\-mean random variables with varianceσ2\\sigma^\{2\}on a sequence of uniform grids over the fixed spatiotemporal domainQQ\. By*grid refinement*we meannt,nx→∞n\_\{t\},n\_\{x\}\\to\\inftyon this fixed domain, withht,hx→0h\_\{t\},h\_\{x\}\\to 0andht=Θ⁡\(nt−1\)h\_\{t\}=\\Theta\(n\_\{t\}^\{\-1\}\),hx=Θ⁡\(nx−1\)h\_\{x\}=\\Theta\(n\_\{x\}^\{\-1\}\)\. Define

⟨f,g⟩h=ht​hx​∑i=0nt−1∑j=0nx−1fi​j​gi​j,∥v∥h2=ht​hx​∑i=0nt−1∑j=0nx−1vi​j2\.\\langle f,\\,g\\rangle\_\{h\}=h\_\{t\}h\_\{x\}\\sum\_\{i=0\}^\{n\_\{t\}\-1\}\\sum\_\{j=0\}^\{n\_\{x\}\-1\}f\_\{ij\}g\_\{ij\},\\hskip 20\.00003pt\\lVert v\\rVert\_\{h\}^\{2\}=h\_\{t\}h\_\{x\}\\sum\_\{i=0\}^\{n\_\{t\}\-1\}\\sum\_\{j=0\}^\{n\_\{x\}\-1\}v\_\{ij\}^\{2\}\.Thus∥⋅∥h\\lVert\\cdot\\rVert\_\{h\}is the quadrature\-weighted discreteL2​\(Q\)L^\{2\}\(Q\)norm on grid samples\. Let∥v∥L2​\(Q\)2=∫Q\|v\|2​𝑑x​𝑑t\\lVert v\\rVert\_\{L^\{2\}\(Q\)\}^\{2\}=\\int\_\{Q\}\|v\|^\{2\}\\,\\,\\mathrm\{d\}x\\,\\mathrm\{d\}t\. For part \(i\), hold the test functionϕk\\phi\_\{k\}fixed onQQas the grid is refined and assume that it belongs to the adjoint domain, so thatωk,βh=𝒳β∗,h​ϕk\\omega^\{h\}\_\{k,\\beta\}=\\mathcal\{X\}\_\{\\beta\}^\{\\ast,h\}\\phi\_\{k\}satisfies∥ωk,βh∥h→∥𝒳β∗​ϕk∥L2​\(Q\)<∞\\lVert\\omega^\{h\}\_\{k,\\beta\}\\rVert\_\{h\}\\to\\lVert\\mathcal\{X\}\_\{\\beta\}^\{\\ast\}\\phi\_\{k\}\\rVert\_\{L^\{2\}\(Q\)\}<\\infty\.

1. \(i\)The weak\-feature perturbation \([13](https://arxiv.org/html/2608.12879#S3.E13)\) satisfies Var⁡\(⟨η,ωk,βh⟩h\)=σ2​ht​hx​∥ωk,βh∥h2=O⁡\(ht​hx\)\.\\operatorname\{Var\}\\bigl\(\\langle\\eta,\\,\\omega^\{h\}\_\{k,\\beta\}\\rangle\_\{h\}\\bigr\)=\\sigma^\{2\}h\_\{t\}h\_\{x\}\\lVert\\omega^\{h\}\_\{k,\\beta\}\\rVert\_\{h\}^\{2\}=O\(h\_\{t\}h\_\{x\}\)\.Hence, at fixed domain size, its variance isO⁡\(\(nt​nx\)−1\)O\(\(n\_\{t\}n\_\{x\}\)^\{\-1\}\)as both grid dimensions are refined\.
2. \(ii\)For a periodic pointwise feature whose multipliersβs\_\{\\beta\}is defined above and satisfies\|sβ​\(κ\)\|=\|κ\|β\|s\_\{\\beta\}\(\\kappa\)\|=\|\\kappa\|^\{\\beta\}, Var⁡\(\(𝒳β​η\)i​j\)=σ2nx​∑ℓ=−⌊nx/2⌋⌈nx/2⌉−1\|sβ​\(κℓ\)\|2∼π2​β2​β\+1​σ2​hx−2​β,κℓ=2​π​ℓLx\.\\operatorname\{Var\}\\bigl\(\(\\mathcal\{X\}\_\{\\beta\}\\eta\)\_\{ij\}\\bigr\)=\\frac\{\\sigma^\{2\}\}\{n\_\{x\}\}\\sum\_\{\\ell=\-\\lfloor n\_\{x\}/2\\rfloor\}^\{\\lceil n\_\{x\}/2\\rceil\-1\}\|s\_\{\\beta\}\(\\kappa\_\{\\ell\}\)\|^\{2\}\\sim\\frac\{\\pi^\{2\\beta\}\}\{2\\beta\+1\}\\sigma^\{2\}h\_\{x\}^\{\-2\\beta\},\\hskip 20\.00003pt\\kappa\_\{\\ell\}=\\frac\{2\\pi\\ell\}\{L\_\{x\}\}\.Thus the variance of every unregularised positive\-order strong\-form feature \(β\>0\\beta\>0\) diverges under spatial refinement, at a rate that increases withβ\\beta; the separately defined identity candidate atβ=0\\beta=0retains varianceσ2\\sigma^\{2\}\. For an even grid, a real\-valued directional implementation may treat the single unpaired Nyquist mode separately; this changes one summand only and leaves the asymptotic relation unchanged, as noted in the proof\.

###### Proof\.

For \(i\), independence removes all cross\-covariances:

Var⁡\(ht​hx​∑i=0nt−1∑j=0nx−1ηi​j​\(ωk,βh\)i​j\)=σ2​\(ht​hx\)2​∑i=0nt−1∑j=0nx−1\(ωk,βh\)i​j2=σ2​ht​hx​∥ωk,βh∥h2\.\\operatorname\{Var\}\\\!\\left\(h\_\{t\}h\_\{x\}\\sum\_\{i=0\}^\{n\_\{t\}\-1\}\\sum\_\{j=0\}^\{n\_\{x\}\-1\}\\eta\_\{ij\}\(\\omega^\{h\}\_\{k,\\beta\}\)\_\{ij\}\\right\)=\\sigma^\{2\}\(h\_\{t\}h\_\{x\}\)^\{2\}\\sum\_\{i=0\}^\{n\_\{t\}\-1\}\\sum\_\{j=0\}^\{n\_\{x\}\-1\}\(\\omega^\{h\}\_\{k,\\beta\}\)\_\{ij\}^\{2\}=\\sigma^\{2\}h\_\{t\}h\_\{x\}\\lVert\\omega^\{h\}\_\{k,\\beta\}\\rVert\_\{h\}^\{2\}\.The assumed norm convergence makes the final factor bounded\. It holds for the fixed smooth Fourier and periodised\-Gaussian tests used with the periodic operators; for a one\-sided operator, it holds whenever the fixed test function belongs to the corresponding adjoint domain\.

For \(ii\), fix a time index and writeηj\\eta\_\{j\}for the spatial noise samples on that time slice\. Let

ℒnx=\{−⌊nx/2⌋,…,⌈nx/2⌉−1\},κℓ=2​π​ℓLx\.\\mathcal\{L\}\_\{n\_\{x\}\}=\\\{\-\\lfloor n\_\{x\}/2\\rfloor,\\ldots,\\lceil n\_\{x\}/2\\rceil\-1\\\},\\hskip 20\.00003pt\\kappa\_\{\\ell\}=\\frac\{2\\pi\\ell\}\{L\_\{x\}\}\.With the discrete Fourier transform

η^ℓ=ℱx​\[η\]​\(κℓ\)=∑q=0nx−1ηq​e−i​κℓ​xq,\\widehat\{\\eta\}\_\{\\ell\}=\\mathcal\{F\}\_\{x\}\[\\eta\]\(\\kappa\_\{\\ell\}\)=\\sum\_\{q=0\}^\{n\_\{x\}\-1\}\\eta\_\{q\}e^\{\-\\mathrm\{i\}\\kappa\_\{\\ell\}x\_\{q\}\},the multiplier definition gives the inverse representation

\(𝒳β​η\)​\(xj\)=1nx​∑ℓ∈ℒnxsβ​\(κℓ\)​η^ℓ​ei​κℓ​xj\.\(\\mathcal\{X\}\_\{\\beta\}\\eta\)\(x\_\{j\}\)=\\frac\{1\}\{n\_\{x\}\}\\sum\_\{\\ell\\in\\mathcal\{L\}\_\{n\_\{x\}\}\}s\_\{\\beta\}\(\\kappa\_\{\\ell\}\)\\widehat\{\\eta\}\_\{\\ell\}e^\{\\mathrm\{i\}\\kappa\_\{\\ell\}x\_\{j\}\}\.Hereδq​r\\delta\_\{qr\}denotes the Kronecker delta\. Since𝔼⁡\[ηq​ηr\]=σ2​δq​r\\mathbb\{E\}\[\\eta\_\{q\}\\eta\_\{r\}\]=\\sigma^\{2\}\\delta\_\{qr\}, discrete Fourier orthogonality yields

𝔼⁡\[η^ℓ​η^m¯\]\\displaystyle\\mathbb\{E\}\[\\widehat\{\\eta\}\_\{\\ell\}\\overline\{\\widehat\{\\eta\}\_\{m\}\}\]=∑q=0nx−1∑r=0nx−1𝔼⁡\[ηq​ηr\]​e−i​κℓ​xq​ei​κm​xr\\displaystyle=\\sum\_\{q=0\}^\{n\_\{x\}\-1\}\\sum\_\{r=0\}^\{n\_\{x\}\-1\}\\mathbb\{E\}\[\\eta\_\{q\}\\eta\_\{r\}\]e^\{\-\\mathrm\{i\}\\kappa\_\{\\ell\}x\_\{q\}\}e^\{\\mathrm\{i\}\\kappa\_\{m\}x\_\{r\}\}=σ2​∑q=0nx−1e−i⁡\(κℓ−κm\)​xq=nx​σ2​δℓ​m\.\\displaystyle=\\sigma^\{2\}\\sum\_\{q=0\}^\{n\_\{x\}\-1\}e^\{\-\\mathrm\{i\}\(\\kappa\_\{\\ell\}\-\\kappa\_\{m\}\)x\_\{q\}\}=n\_\{x\}\\sigma^\{2\}\\delta\_\{\\ell m\}\.The output is real for the conjugate\-symmetric multipliers considered here\. Its mean is zero, and expanding the squared magnitude therefore gives

Var⁡\(\(𝒳β​η\)​\(xj\)\)\\displaystyle\\operatorname\{Var\}\\bigl\(\(\\mathcal\{X\}\_\{\\beta\}\\eta\)\(x\_\{j\}\)\\bigr\)=1nx2​∑ℓ,m∈ℒnxsβ​\(κℓ\)​sβ​\(κm\)¯​ei⁡\(κℓ−κm\)​xj​𝔼​\[η^ℓ​η^m¯\]\\displaystyle=\\frac\{1\}\{n\_\{x\}^\{2\}\}\\sum\_\{\\ell,m\\in\\mathcal\{L\}\_\{n\_\{x\}\}\}s\_\{\\beta\}\(\\kappa\_\{\\ell\}\)\\overline\{s\_\{\\beta\}\(\\kappa\_\{m\}\)\}e^\{\\mathrm\{i\}\(\\kappa\_\{\\ell\}\-\\kappa\_\{m\}\)x\_\{j\}\}\\mathbb\{E\}\[\\widehat\{\\eta\}\_\{\\ell\}\\overline\{\\widehat\{\\eta\}\_\{m\}\}\]=σ2nx​∑ℓ=−⌊nx/2⌋⌈nx/2⌉−1\|sβ​\(κℓ\)\|2=σ2nx​∑ℓ=−⌊nx/2⌋⌈nx/2⌉−1\|κℓ\|2​β\.\\displaystyle=\\frac\{\\sigma^\{2\}\}\{n\_\{x\}\}\\sum\_\{\\ell=\-\\lfloor n\_\{x\}/2\\rfloor\}^\{\\lceil n\_\{x\}/2\\rceil\-1\}\|s\_\{\\beta\}\(\\kappa\_\{\\ell\}\)\|^\{2\}=\\frac\{\\sigma^\{2\}\}\{n\_\{x\}\}\\sum\_\{\\ell=\-\\lfloor n\_\{x\}/2\\rfloor\}^\{\\lceil n\_\{x\}/2\\rceil\-1\}\|\\kappa\_\{\\ell\}\|^\{2\\beta\}\.Becausehx=Lx/nxh\_\{x\}=L\_\{x\}/n\_\{x\}, the last sum obeys

hx2​βnx​∑ℓ=−⌊nx/2⌋⌈nx/2⌉−1\|κℓ\|2​β\\displaystyle\\frac\{h\_\{x\}^\{2\\beta\}\}\{n\_\{x\}\}\\sum\_\{\\ell=\-\\lfloor n\_\{x\}/2\\rfloor\}^\{\\lceil n\_\{x\}/2\\rceil\-1\}\|\\kappa\_\{\\ell\}\|^\{2\\beta\}=1nx​∑ℓ=−⌊nx/2⌋⌈nx/2⌉−1\|2​π​ℓnx\|2​β\\displaystyle=\\frac\{1\}\{n\_\{x\}\}\\sum\_\{\\ell=\-\\lfloor n\_\{x\}/2\\rfloor\}^\{\\lceil n\_\{x\}/2\\rceil\-1\}\\left\|\\frac\{2\\pi\\ell\}\{n\_\{x\}\}\\right\|^\{2\\beta\}⟶∫−1/21/2\|2πr\|2​βdr=π2​β2​β\+1\.\\displaystyle\\longrightarrow\\int\_\{\-1/2\}^\{1/2\}\|2\\pi r\|^\{2\\beta\}\\,\\,\\mathrm\{d\}r=\\frac\{\\pi^\{2\\beta\}\}\{2\\beta\+1\}\.Consequently,

Var⁡\(\(𝒳β​η\)​\(xj\)\)∼π2​β2​β\+1​σ2​hx−2​β,\\operatorname\{Var\}\\bigl\(\(\\mathcal\{X\}\_\{\\beta\}\\eta\)\(x\_\{j\}\)\\bigr\)\\sim\\frac\{\\pi^\{2\\beta\}\}\{2\\beta\+1\}\\sigma^\{2\}h\_\{x\}^\{\-2\\beta\},which is exactly the asymptotic relation stated in part \(ii\)\. For the even\-grid directional implementation, the real\-valued projection affects only the unpaired Nyquist mode\. Its contribution to the variance isO⁡\(hx−\(2​β−1\)\)O\(h\_\{x\}^\{\-\(2\\beta\-1\)\}\), one order lower than theO⁡\(hx−2​β\)O\(h\_\{x\}^\{\-2\\beta\}\)total, and therefore does not change the leading constant or rate\. ∎

Proposition[1](https://arxiv.org/html/2608.12879#Thmtheorem1)gives a feature\-level explanation of why denser observations benefit the linear weak formulation more than the unregularised strong form under its assumptions\. With the implementation’s separate discreteℓ2\\ell^\{2\}row normalisation, the variance isO⁡\(ht2​hx2\)O\(h\_\{t\}^\{2\}h\_\{x\}^\{2\}\)for fixed separable tests, while signal and noise are rescaled identically\. On a fixed domain, each linear weak measurement averages more independent noise samples and becomes more stable, whereas pointwise fractional differentiation admits progressively higher wavenumbers and amplifies their noise\. Nonetheless, correlated noise, the Caputo initial\-data terms, nonlinear features, and conditioning of the regression can alter the final estimator\. We therefore assess coefficient and structure recovery empirically in Section[5](https://arxiv.org/html/2608.12879#S5), and Appendix[8](https://arxiv.org/html/2608.12879#S8)verifies the predicted rates numerically\.

###### Corollary 2\(Multiplicative measurement noise\)\.

Let the measured field beui​j=ui​j⋆​\(1\+ρ​ζi​j\)u\_\{ij\}=u^\{\\star\}\_\{ij\}\(1\+\\rho\\zeta\_\{ij\}\), where theζi​j\\zeta\_\{ij\}are independent, have zero mean and varianceσζ2\\sigma\_\{\\zeta\}^\{2\}, andu⋆u^\{\\star\}is regarded as fixed\. For the same fixed test function onQQand discrete weak\-feature weightsωk,βh\\omega^\{h\}\_\{k,\\beta\}as in Proposition[1](https://arxiv.org/html/2608.12879#Thmtheorem1)\(i\), the conditional perturbation variance is

Var⁡\(⟨u−u⋆,ωk,βh⟩h\|u⋆\)\\displaystyle\\operatorname\{Var\}\\\!\\left\(\\langle u\-u^\{\\star\},\\,\\omega^\{h\}\_\{k,\\beta\}\\rangle\_\{h\}\\,\\middle\|\\,u^\{\\star\}\\right\)=ρ2​σζ2​\(ht​hx\)2​∑i=0nt−1∑j=0nx−1\(ui​j⋆​\(ωk,βh\)i​j\)2\\displaystyle=\\rho^\{2\}\\sigma\_\{\\zeta\}^\{2\}\(h\_\{t\}h\_\{x\}\)^\{2\}\\sum\_\{i=0\}^\{n\_\{t\}\-1\}\\sum\_\{j=0\}^\{n\_\{x\}\-1\}\\bigl\(u^\{\\star\}\_\{ij\}\(\\omega^\{h\}\_\{k,\\beta\}\)\_\{ij\}\\bigr\)^\{2\}≤ρ2​σζ2​∥u⋆∥∞2​ht​hx​∥ωk,βh∥h2=O⁡\(ht​hx\)\.\\displaystyle\\leq\\rho^\{2\}\\sigma\_\{\\zeta\}^\{2\}\\lVert u^\{\\star\}\\rVert\_\{\\infty\}^\{2\}h\_\{t\}h\_\{x\}\\lVert\\omega^\{h\}\_\{k,\\beta\}\\rVert\_\{h\}^\{2\}=O\(h\_\{t\}h\_\{x\}\)\.For the uniform perturbationsζi​j∼𝒰⁡\[−1,1\]\\zeta\_\{ij\}\\sim\\mathcal\{U\}\[\-1,1\]used in the main experiments,σζ2=1/3\\sigma\_\{\\zeta\}^\{2\}=1/3\.

###### Proof\.

Condition on the noise\-free fieldu⋆u^\{\\star\}\. The perturbation is

⟨u−u⋆,ωk,βh⟩h=ρ​ht​hx​∑i=0nt−1∑j=0nx−1ui​j⋆​ζi​j​\(ωk,βh\)i​j\.\\langle u\-u^\{\\star\},\\,\\omega^\{h\}\_\{k,\\beta\}\\rangle\_\{h\}=\\rho h\_\{t\}h\_\{x\}\\sum\_\{i=0\}^\{n\_\{t\}\-1\}\\sum\_\{j=0\}^\{n\_\{x\}\-1\}u^\{\\star\}\_\{ij\}\\zeta\_\{ij\}\(\\omega^\{h\}\_\{k,\\beta\}\)\_\{ij\}\.Its conditional mean is zero\. Independence again removes the cross\-terms, giving

Var⁡\(⟨u−u⋆,ωk,βh⟩h\|u⋆\)=ρ2​σζ2​\(ht​hx\)2​∑i=0nt−1∑j=0nx−1\(ui​j⋆\)2​\(ωk,βh\)i​j2\.\\operatorname\{Var\}\\\!\\left\(\\langle u\-u^\{\\star\},\\,\\omega^\{h\}\_\{k,\\beta\}\\rangle\_\{h\}\\,\\middle\|\\,u^\{\\star\}\\right\)=\\rho^\{2\}\\sigma\_\{\\zeta\}^\{2\}\(h\_\{t\}h\_\{x\}\)^\{2\}\\sum\_\{i=0\}^\{n\_\{t\}\-1\}\\sum\_\{j=0\}^\{n\_\{x\}\-1\}\(u^\{\\star\}\_\{ij\}\)^\{2\}\(\\omega^\{h\}\_\{k,\\beta\}\)\_\{ij\}^\{2\}\.Using\|ui​j⋆\|≤∥u⋆∥∞\|u^\{\\star\}\_\{ij\}\|\\leq\\lVert u^\{\\star\}\\rVert\_\{\\infty\}and the definition of∥⋅∥h\\lVert\\cdot\\rVert\_\{h\}gives, explicitly,

Var⁡\(⟨u−u⋆,ωk,βh⟩h\|u⋆\)\\displaystyle\\operatorname\{Var\}\\\!\\left\(\\langle u\-u^\{\\star\},\\,\\omega^\{h\}\_\{k,\\beta\}\\rangle\_\{h\}\\,\\middle\|\\,u^\{\\star\}\\right\)≤ρ2​σζ2​∥u⋆∥∞2​\(ht​hx\)2​∑i=0nt−1∑j=0nx−1\(ωk,βh\)i​j2\\displaystyle\\leq\\rho^\{2\}\\sigma\_\{\\zeta\}^\{2\}\\lVert u^\{\\star\}\\rVert\_\{\\infty\}^\{2\}\(h\_\{t\}h\_\{x\}\)^\{2\}\\sum\_\{i=0\}^\{n\_\{t\}\-1\}\\sum\_\{j=0\}^\{n\_\{x\}\-1\}\(\\omega^\{h\}\_\{k,\\beta\}\)\_\{ij\}^\{2\}=ρ2​σζ2​∥u⋆∥∞2​ht​hx​∥ωk,βh∥h2\\displaystyle=\\rho^\{2\}\\sigma\_\{\\zeta\}^\{2\}\\lVert u^\{\\star\}\\rVert\_\{\\infty\}^\{2\}h\_\{t\}h\_\{x\}\\lVert\\omega^\{h\}\_\{k,\\beta\}\\rVert\_\{h\}^\{2\}=O⁡\(ht​hx\),\\displaystyle=O\(h\_\{t\}h\_\{x\}\),where the final step uses the bounded\-weight assumption from Proposition[1](https://arxiv.org/html/2608.12879#Thmtheorem1)\(i\)\. ∎

##### Subunit branch\.

The same noisy initial valueu⁡\(0,xj\)u\(0,x\_\{j\}\)is reused in every temporal contribution\. Its variance term is

σ2​ht2​hx2​∑j=0nx−1\(∑i=0nt−1ωi​j\)2\.\\sigma^\{2\}h\_\{t\}^\{2\}h\_\{x\}^\{2\}\\sum\_\{j=0\}^\{n\_\{x\}\-1\}\\left\(\\sum\_\{i=0\}^\{n\_\{t\}\-1\}\\omega\_\{ij\}\\right\)^\{2\}\.For samples of a fixed continuum test weight, an ordinary full\-field contribution has∑i,jωi​j2=O⁡\(\(ht​hx\)−1\)\\sum\_\{i,j\}\\omega\_\{ij\}^\{2\}=O\(\(h\_\{t\}h\_\{x\}\)^\{\-1\}\), and multiplication byht2​hx2h\_\{t\}^\{2\}h\_\{x\}^\{2\}gives theO⁡\(ht​hx\)O\(h\_\{t\}h\_\{x\}\)variance in Proposition[1](https://arxiv.org/html/2608.12879#Thmtheorem1)\. Reusing the initial sample instead gives

∑j=0nx−1\(∑i=0nt−1ωi​j\)2=O⁡\(ht−2​hx−1\),\\sum\_\{j=0\}^\{n\_\{x\}\-1\}\\left\(\\sum\_\{i=0\}^\{n\_\{t\}\-1\}\\omega\_\{ij\}\\right\)^\{2\}=O\(h\_\{t\}^\{\-2\}h\_\{x\}^\{\-1\}\),and hence a variance ofO⁡\(hx\)O\(h\_\{x\}\)\. The squared quadrature factors are therefore offset by the growth of the discrete sums\. Spatial refinement still averages independent initial\-time samples, but temporal refinement does not create new independent copies of the initial datum\. This slower rate does not by itself determine the direction of an order\-selection bias\.

##### Superunit branch\.

For1<α<21<\\alpha<2, the reported evaluator usesDhα=𝖫𝟣α−1​D1D\_\{h\}^\{\\alpha\}=\\mathsf\{L1\}\_\{\\alpha\-1\}D\_\{1\}\. For the separable rowϕk​\(t,x\)=ϑℓt​\(t\)​ψℓx​\(x\)\\phi\_\{k\}\(t,x\)=\\vartheta\_\{\\ell\_\{t\}\}\(t\)\\psi\_\{\\ell\_\{x\}\}\(x\), consider unnormalised samples of the fixed continuum tests and define

𝒘\(α\)=\(Dhα\)⊤ϑℓt,ωi​j=wi\(α\)\(𝝍ℓx\)j,i=0,…,nt−1,j=0,…,nx−1\.\\bm\{w\}^\{\(\\alpha\)\}=\(D\_\{h\}^\{\\alpha\}\)^\{\\top\}\\bm\{\\vartheta\}\_\{\\ell\_\{t\}\},\\hskip 20\.00003pt\\omega\_\{ij\}=w\_\{i\}^\{\(\\alpha\)\}\(\\bm\{\\psi\}\_\{\\ell\_\{x\}\}\)\_\{j\},\\hskip 10\.00002pti=0,\\ldots,n\_\{t\}\-1,\\hskip 10\.00002ptj=0,\\ldots,n\_\{x\}\-1\.The factorisation also explains the endpoint weights\. Set𝒗=𝖫𝟣α−1⊤​ϑℓt\\bm\{v\}=\\mathsf\{L1\}\_\{\\alpha\-1\}^\{\\top\}\\bm\{\\vartheta\}\_\{\\ell\_\{t\}\}; then𝒘\(α\)=D1⊤​𝒗\\bm\{w\}^\{\(\\alpha\)\}=D\_\{1\}^\{\\top\}\\bm\{v\}\. For the fixed smooth temporal tests used under refinement, the L1 weights givev0=O⁡\(ht−1\)v\_\{0\}=O\(h\_\{t\}^\{\-1\}\), while the neighbouring entriesv1v\_\{1\}andv2v\_\{2\}remainO⁡\(1\)O\(1\)\. Because the first row ofD1D\_\{1\}is a forward difference and the next rows use centred differences, the first two transpose weights satisfyw0\(α\)=−v0/ht−v1/\(2ht\)w\_\{0\}^\{\(\\alpha\)\}=\-v\_\{0\}/h\_\{t\}\-v\_\{1\}/\(2h\_\{t\}\)andw1\(α\)=v0/ht−v2/\(2​ht\)w\_\{1\}^\{\(\\alpha\)\}=v\_\{0\}/h\_\{t\}\-v\_\{2\}/\(2h\_\{t\}\)\. Hencew0\(α\),w1\(α\)=O⁡\(ht−2\)w\_\{0\}^\{\(\\alpha\)\},w\_\{1\}^\{\(\\alpha\)\}=O\(h\_\{t\}^\{\-2\}\)andw0\(α\)=−w1\(α\)\+O⁡\(ht−1\)w\_\{0\}^\{\(\\alpha\)\}=\-w\_\{1\}^\{\(\\alpha\)\}\+O\(h\_\{t\}^\{\-1\}\)\. This opposite\-sign initial pair therefore dominates∑i\|wi\(α\)\|2=O⁡\(ht−4\)\\sum\_\{i\}\|w\_\{i\}^\{\(\\alpha\)\}\|^\{2\}=O\(h\_\{t\}^\{\-4\}\)\. Moreover,

hx​∑j=0nx−1\|\(𝝍ℓx\)j\|2⟶∥ψℓx∥L2​\(Ω\)2,h\_\{x\}\\sum\_\{j=0\}^\{n\_\{x\}\-1\}\|\(\\bm\{\\psi\}\_\{\\ell\_\{x\}\}\)\_\{j\}\|^\{2\}\\longrightarrow\\lVert\\psi\_\{\\ell\_\{x\}\}\\rVert\_\{L^\{2\}\(\\Omega\)\}^\{2\},hence∑j\|\(𝝍ℓx\)j\|2=O⁡\(hx−1\)\\sum\_\{j\}\|\(\\bm\{\\psi\}\_\{\\ell\_\{x\}\}\)\_\{j\}\|^\{2\}=O\(h\_\{x\}^\{\-1\}\)\. For independent additive noise, the unnormalised target\-row variance is therefore

σ2​ht2​hx2​\(∑i=0nt−1\|wi\(α\)\|2\)​\(∑j=0nx−1\|\(𝝍ℓx\)j\|2\)=O⁡\(σ2​hx​ht−2\)\.\\sigma^\{2\}h\_\{t\}^\{2\}h\_\{x\}^\{2\}\\left\(\\sum\_\{i=0\}^\{n\_\{t\}\-1\}\|w\_\{i\}^\{\(\\alpha\)\}\|^\{2\}\\right\)\\left\(\\sum\_\{j=0\}^\{n\_\{x\}\-1\}\|\(\\bm\{\\psi\}\_\{\\ell\_\{x\}\}\)\_\{j\}\|^\{2\}\\right\)=O\\\!\\left\(\\sigma^\{2\}h\_\{x\}h\_\{t\}^\{\-2\}\\right\)\.For the multiplicative model of Corollary[2](https://arxiv.org/html/2608.12879#Thmtheorem2)with boundedu⋆u^\{\\star\}, the same argument gives the upper boundO⁡\(ρ2​σζ2​∥u⋆∥∞2​hx​ht−2\)O\(\\rho^\{2\}\\sigma\_\{\\zeta\}^\{2\}\\lVert u^\{\\star\}\\rVert\_\{\\infty\}^\{2\}h\_\{x\}h\_\{t\}^\{\-2\}\)\. Thus the endpoint sensitivity also applies to the noise law used in the superunit diagnostic\. Under the implementation’s separate discreteℓ2\\ell^\{2\}normalisation of the temporal and spatial test rows, the corresponding absolute variance isO⁡\(σ2​hx2​ht−1\)O\(\\sigma^\{2\}h\_\{x\}^\{2\}h\_\{t\}^\{\-1\}\)\. The clean target is rescaled by the same factors; this normalisation therefore does not remove the relative sensitivity to initial\-time noise\.

The practical implication is that the discrete transpose removes pointwise interior differentiation while retaining the initial\-boundary dependence intrinsic to the Caputo derivative\. Ordinary endpoint vanishing ofϑℓt\\vartheta\_\{\\ell\_\{t\}\}leaves the initial pair coupled through𝖫𝟣α−1⊤\\mathsf\{L1\}\_\{\\alpha\-1\}^\{\\top\}and the first rows ofD1⊤D\_\{1\}^\{\\top\}\. This behaviour is consistent with the initial\-rate trace in the continuous identity\. The complete target also contains the noisy full\-field projection and is scored after variance normalisation in Eq\. \([21](https://arxiv.org/html/2608.12879#S4.E21)\)\.

Appendix[8](https://arxiv.org/html/2608.12879#S8)separates these effects with five fixed\-active\-set cases: fully clean data; noise in every data path; a noisy field with the initial slice restored to its clean value; noise only in the temporal target; and noise only in the right\-hand\-side library\. Comparing these cases identifies which data path drives a selected\-order shift\.

Corollary[2](https://arxiv.org/html/2608.12879#Thmtheorem2)covers the principal experimental noise model: although multiplicative noise is heteroscedastic, the variance of a fixed linear weak feature still vanishes under refinement\. Appendix[13](https://arxiv.org/html/2608.12879#S13)shows the same weak–strong separation under additive Gaussian noise\.

The weak formulation reduces noise amplification, but accurate order recovery also requires the fractional order to be distinguishable within the weak regression\. The following subsection quantifies this second issue for spatial fractional orders\.

#### 3\.5Local identifiability of spatial fractional orders

On a fixed active support, we call a spatial orderβj\\beta\_\{j\}locally identifiable when a small change inβj\\beta\_\{j\}produces a change in the weak regression that cannot be reproduced by refitting the active coefficients\. This distinction matters because a weak formulation can suppress measurement noise while nearby orders still generate very similar regression columns\. To isolate the order effect, fix the temporal branch, the active support, and all continuous orders exceptβj\\beta\_\{j\}\. Using exact \(non\-interpolated\) weak features, write

𝒃=Θ⁡\(βj\)​𝝃\+𝒓,Θ=\[𝜽1,…,𝜽c\]∈ℝK×c\.\\bm\{b\}=\\Theta\(\\beta\_\{j\}\)\\bm\{\\xi\}\+\\bm\{r\},\\hskip 20\.00003pt\\Theta=\[\\bm\{\\theta\}\_\{1\},\\ldots,\\bm\{\\theta\}\_\{c\}\]\\in\\mathbb\{R\}^\{K\\times c\}\.\(14\)At the reference order,𝜽˙j=∂𝜽j/∂βj\\dot\{\\bm\{\\theta\}\}\_\{j\}=\\partial\\bm\{\\theta\}\_\{j\}/\\partial\\beta\_\{j\}measures how thejjth weak column changes withβj\\beta\_\{j\}\. LetPΘ=Θ​Θ†P\_\{\\Theta\}=\\Theta\\Theta^\{\\dagger\}denote the orthogonal projector ontocol⁡\(Θ\)\\operatorname\{col\}\(\\Theta\), whereΘ†\\Theta^\{\\dagger\}is the Moore–Penrose pseudoinverse\.

###### Proposition 3\(Local order sensitivity after coefficient refitting\)\.

Suppose𝛉j≠0\\bm\{\\theta\}\_\{j\}\\neq 0andξj≠0\\xi\_\{j\}\\neq 0\. For an infinitesimal perturbationβj↦βj\+δ​βj\\beta\_\{j\}\\mapsto\\beta\_\{j\}\+\\delta\\beta\_\{j\}, allow the coefficients to refit as𝛏↦𝛏\+δ​βj​𝐯\\bm\{\\xi\}\\mapsto\\bm\{\\xi\}\+\\delta\\beta\_\{j\}\\bm\{v\}\. Then the smallest first\-order change in the model mean is

min𝒗∈ℝc⁡‖Θ​𝒗\+ξj​𝜽˙j‖2=\|ξj\|​‖\(I−PΘ\)​𝜽˙j‖2\.\\min\_\{\\bm\{v\}\\in\\mathbb\{R\}^\{c\}\}\\left\\lVert\\Theta\\bm\{v\}\+\\xi\_\{j\}\\dot\{\\bm\{\\theta\}\}\_\{j\}\\right\\rVert\_\{2\}=\|\\xi\_\{j\}\|\\,\\left\\lVert\(I\-P\_\{\\Theta\}\)\\dot\{\\bm\{\\theta\}\}\_\{j\}\\right\\rVert\_\{2\}\.\(15\)Henceβj\\beta\_\{j\}is first\-order confounded with coefficient changes when𝛉˙j∈col⁡\(Θ\)\\dot\{\\bm\{\\theta\}\}\_\{j\}\\in\\operatorname\{col\}\(\\Theta\)\.

###### Proof\.

Taylor expansion gives

Θ⁡\(βj\+δ​βj\)​\(𝝃\+δ​βj​𝒗\)−Θ⁡\(βj\)​𝝃=δ​βj​\(Θ​𝒗\+ξj​𝜽˙j\)\+o⁡\(\|δ​βj\|\)\.\\Theta\(\\beta\_\{j\}\+\\delta\\beta\_\{j\}\)\(\\bm\{\\xi\}\+\\delta\\beta\_\{j\}\\bm\{v\}\)\-\\Theta\(\\beta\_\{j\}\)\\bm\{\\xi\}=\\delta\\beta\_\{j\}\\bigl\(\\Theta\\bm\{v\}\+\\xi\_\{j\}\\dot\{\\bm\{\\theta\}\}\_\{j\}\\bigr\)\+o\(\|\\delta\\beta\_\{j\}\|\)\.The termξj​𝜽˙j\\xi\_\{j\}\\dot\{\\bm\{\\theta\}\}\_\{j\}is the first\-order change caused by perturbing the order\. Refitting the coefficients contributesΘ​𝒗\\Theta\\bm\{v\}, which can reproduce any vector incol⁡\(Θ\)\\operatorname\{col\}\(\\Theta\)\. Least squares therefore cancels the component ofξj​𝜽˙j\\xi\_\{j\}\\dot\{\\bm\{\\theta\}\}\_\{j\}that lies in this column space\. The component that remains isξj​\(I−PΘ\)​𝜽˙j\\xi\_\{j\}\(I\-P\_\{\\Theta\}\)\\dot\{\\bm\{\\theta\}\}\_\{j\}; taking its norm gives Eq\. \([15](https://arxiv.org/html/2608.12879#S3.E15)\)\. ∎

This result motivates the dimensionless, coefficient\-independent diagnostic

Sβj=‖\(I−PΘ\)​𝜽˙j‖2∥𝜽j∥2\.S\_\{\\beta\_\{j\}\}=\\frac\{\\left\\lVert\(I\-P\_\{\\Theta\}\)\\dot\{\\bm\{\\theta\}\}\_\{j\}\\right\\rVert\_\{2\}\}\{\\lVert\\bm\{\\theta\}\_\{j\}\\rVert\_\{2\}\}\.\(16\)The numerator measures the part of the column change that remains after the best first\-order coefficient adjustment, and the denominator removes the scale of the column itself\. ThusSβj=0S\_\{\\beta\_\{j\}\}=0corresponds to complete first\-order confounding on the fixed support, whereas a larger value indicates that the order produces a more distinct regression direction\. The diagnostic is local and support\-conditioned; in Section[5\.2](https://arxiv.org/html/2608.12879#S5.SS2)we compare it with fixed\-support order errors under noise\.

A simple Fourier calculation explains why this sensitivity depends on both the operator and the spectral content of the observed field\. Consider one isolated periodic linear term before weak projection, and let

E⁡\(κ\)=∑i=0nt−1\|u^​\(ti,κ\)\|2E\(\\kappa\)=\\sum\_\{i=0\}^\{n\_\{t\}\-1\}\|\\widehat\{u\}\(t\_\{i\},\\kappa\)\|^\{2\}be the total energy carried by spatial wavenumberκ\\kappaover the observed times\. This quantity is relevant because an order can only be distinguished through wavenumbers that are actually present in the data\. Both periodic operators satisfy\|sβ​\(κ\)\|=\|κ\|β\|s\_\{\\beta\}\(\\kappa\)\|=\|\\kappa\|^\{\\beta\}, so the fraction of operator\-weighted energy carried by each nonzero wavenumber is

qβ​\(κ\)=\|κ\|2​β​E​\(κ\)∑ν≠0\|ν\|2​β​E​\(ν\),κ≠0\.q\_\{\\beta\}\(\\kappa\)=\\frac\{\|\\kappa\|^\{2\\beta\}E\(\\kappa\)\}\{\\sum\_\{\\nu\\neq 0\}\|\\nu\|^\{2\\beta\}E\(\\nu\)\},\\hskip 20\.00003pt\\kappa\\neq 0\.Applying the same coefficient\-refitting argument to the Fourier multipliers gives the isolated\-term sensitivity

Sβ,spec2=\{Varqβ⁡\[log⁡\|κ\|\],Riesz,Varqβ⁡\[log⁡\|κ\|\]\+π2/4,directional\.S\_\{\\beta,\\mathrm\{spec\}\}^\{2\}=\\begin\{cases\}\\operatorname\{Var\}\_\{q\_\{\\beta\}\}\[\\log\|\\kappa\|\],&\\text\{Riesz\},\\\\\[2\.0pt\] \\operatorname\{Var\}\_\{q\_\{\\beta\}\}\[\\log\|\\kappa\|\]\+\\pi^\{2\}/4,&\\text\{directional\}\.\\end\{cases\}\(17\)HereVarqβ⁡\[log⁡\|κ\|\]\\operatorname\{Var\}\_\{q\_\{\\beta\}\}\[\\log\|\\kappa\|\]measures the spread of the operator\-weighted wavenumbers on a logarithmic scale\. For a Riesz term, this spread is the entire source of local order sensitivity: if all spectral weight lies at one\|κ\|\|\\kappa\|, changingβ\\betaonly rescales the feature and can be absorbed exactly by the coefficient\. With energy at several distinct wavenumbers, changingβ\\betaalters their relative magnitudes and becomes easier to distinguish\. A directional derivative has the same magnitude scaling and also changes phase throughexp⁡\(i​π​β​sgn⁡\(κ\)/2\)\\exp\(\\mathrm\{i\}\\pi\\beta\\operatorname\{sgn\}\(\\kappa\)/2\), which contributes the additionalπ2/4\\pi^\{2\}/4term for a real field with conjugate\-symmetric spectral energy\. Eq\. \([17](https://arxiv.org/html/2608.12879#S3.E17)\) therefore explains the operator\-level mechanism\. The reported diagnostic usesSβjS\_\{\\beta\_\{j\}\}in Eq\. \([16](https://arxiv.org/html/2608.12879#S3.E16)\) on the full weak design, so it also incorporates the test functions and the other active columns\.

### 4Pareto\-based subset selection

The second core contribution is a search procedure that treats term type and fractional order jointly without constructing a dense fixed dictionary\.

Algorithm[1](https://arxiv.org/html/2608.12879#alg1)summarises the mechanism of Weak\-Pareto from building the weak feature library to the exact\-order refit\. This section details its search and selection stages\.

#### 4\.1The weak regression system

Stacking theKKweak equations \([7](https://arxiv.org/html/2608.12879#S3.E7)\), one row per test function, gives the following system of equations:

𝒃⁡\(mα,α\)=Θ⁡\(𝒑,𝜷\)​𝝃\+𝒓,bk​\(mα,α\)=⟨𝒯mα,α​u,ϕk⟩,Θk,j=⟨upj​𝒳βj​u,ϕk⟩\.\\bm\{b\}\(m\_\{\\alpha\},\\alpha\)=\\Theta\(\\bm\{p\},\\bm\{\\beta\}\)\\,\\bm\{\\xi\}\+\\bm\{r\},\\hskip 20\.00003ptb\_\{k\}\(m\_\{\\alpha\},\\alpha\)=\\langle\\mathcal\{T\}\_\{m\_\{\\alpha\},\\alpha\}u,\\,\\phi\_\{k\}\\rangle,\\hskip 10\.00002pt\\Theta\_\{k,j\}=\\langle u^\{p\_\{j\}\}\\mathcal\{X\}\_\{\\beta\_\{j\}\}u,\\,\\phi\_\{k\}\\rangle\.\(18\)The target𝒃\\bm\{b\}is assembled from Eq\. \([9](https://arxiv.org/html/2608.12879#S3.E9)\)\. Each design column is constructed from Eq\. \([10](https://arxiv.org/html/2608.12879#S3.E10)\) or \([11](https://arxiv.org/html/2608.12879#S3.E11)\)\. For fixed powers and orders, Eq\. \([18](https://arxiv.org/html/2608.12879#S4.E18)\) is linear in𝝃\\bm\{\\xi\}\. Weak\-Pareto therefore solves the coefficients directly and reserves global optimisation for the fractional orders\.

#### 4\.2Best\-subset regression

Section[2\.2](https://arxiv.org/html/2608.12879#S2.SS2)represents a candidate asℳ=\(mα,α,𝒑,𝜷,𝝃\)\\mathcal\{M\}=\(m\_\{\\alpha\},\\alpha,\\bm\{p\},\\bm\{\\beta\},\\bm\{\\xi\}\)\. At a fixed support sizecc, choosing the discrete temporal branchmαm\_\{\\alpha\}and power vector𝒑\\bm\{p\}leaves the temporal and spatial orders as continuous variables, while the coefficients are obtained by linear regression\. Because the encoding does not predefine an order dictionary, discovery therefore becomes a sequence of best\-subset problems indexed bycc\. Let𝒜=\{0,…,P\}\\mathcal\{A\}=\\\{0,\\dots,P\\\}be the admitted powers\. Itscc\-fold Cartesian product is𝒜c=𝒜×⋯×𝒜\\mathcal\{A\}^\{c\}=\\mathcal\{A\}\\times\\cdots\\times\\mathcal\{A\}, and we retain only the nondecreasing power patterns

𝔓c=\{\(p1,…,pc\)∈𝒜c:p1≤⋯≤pc\}\.\\mathfrak\{P\}\_\{c\}=\\\{\(p\_\{1\},\\dots,p\_\{c\}\)\\in\\mathcal\{A\}^\{c\}:p\_\{1\}\\leq\\cdots\\leq p\_\{c\}\\\}\.For example, if𝒜=\{0,1,2\}\\mathcal\{A\}=\\\{0,1,2\\\}andc=2c=2, then𝔓2=\{\(0,0\),\(0,1\),\(0,2\),\(1,1\),\(1,2\),\(2,2\)\}\\mathfrak\{P\}\_\{2\}=\\\{\(0,0\),\(0,1\),\(0,2\),\(1,1\),\(1,2\),\(2,2\)\\\}\. The ordering removes duplicate permutations of the powers; the associated spatial orders remain continuous and are optimised independently for the terms\. For the declared temporal range\[αmin,αmax\]\[\\alpha\_\{\\min\},\\alpha\_\{\\max\}\]and separatorϵα\\epsilon\_\{\\alpha\}, define

ℐsub\\displaystyle\\mathcal\{I\}\_\{\\mathrm\{sub\}\}=\[αmin,αmax\]∩\(0,1−ϵα\],\\displaystyle=\[\\alpha\_\{\\min\},\\alpha\_\{\\max\}\]\\cap\(0,1\-\\epsilon\_\{\\alpha\}\],ℐint\\displaystyle\\mathcal\{I\}\_\{\\mathrm\{int\}\}=\[αmin,αmax\]∩\{1\},\\displaystyle=\[\\alpha\_\{\\min\},\\alpha\_\{\\max\}\]\\cap\\\{1\\\},ℐsup\\displaystyle\\mathcal\{I\}\_\{\\mathrm\{sup\}\}=\[αmin,αmax\]∩\[1\+ϵα,2\)\.\\displaystyle=\[\\alpha\_\{\\min\},\\alpha\_\{\\max\}\]\\cap\[1\+\\epsilon\_\{\\alpha\},2\)\.The sets correspond directly to the temporal label in the encoding:ℐsub\\mathcal\{I\}\_\{\\mathrm\{sub\}\}contains admitted fractional orders below one,ℐint\\mathcal\{I\}\_\{\\mathrm\{int\}\}contains the exact integer candidateα=1\\alpha=1, andℐsup\\mathcal\{I\}\_\{\\mathrm\{sup\}\}contains admitted fractional orders above one; empty sets are omitted\. For eachc=1,…,cmaxc=1,\\dots,c\_\{\\max\}, we solve

ℳc⋆=arg​min𝐩∈𝔓c,mα,α∈ℐmα,𝜷⁡Jc​\(𝐩,mα,α,𝜷\),Jc=log10⁡\(ℰval\+ε\)\+Πdup,\\mathcal\{M\}\_\{c\}^\{\\star\}=\\argmin\_\{\\bm\{p\}\\in\\mathfrak\{P\}\_\{c\},\\,m\_\{\\alpha\},\\,\\alpha\\in\\mathcal\{I\}\_\{m\_\{\\alpha\}\},\\,\\bm\{\\beta\}\}J\_\{c\}\(\\bm\{p\},m\_\{\\alpha\},\\alpha,\\bm\{\\beta\}\),\\hskip 20\.00003ptJ\_\{c\}=\\log\_\{10\}\\\!\\bigl\(\\mathcal\{E\}\_\{\\mathrm\{val\}\}\+\\varepsilon\\bigr\)\+\\Pi\_\{\\mathrm\{dup\}\},\(19\)whereℰval\\mathcal\{E\}\_\{\\mathrm\{val\}\}is the variance\-normalised held\-out score in Eq\. \([21](https://arxiv.org/html/2608.12879#S4.E21)\) below, andΠdup\\Pi\_\{\\mathrm\{dup\}\}penalises nearly duplicate terms with the same power\. Each branch is optimised independently and scored by the same normalised criterion; changes in the scale of the branch\-specific target therefore do not distort the comparison\. Coefficients are fitted only on the training weak rows\. WithS=diag⁡\(∥𝜽1,tr∥2,…,∥𝜽c,tr∥2\)S=\\operatorname\{diag\}\(\\lVert\\bm\{\\theta\}\_\{1,\\mathrm\{tr\}\}\\rVert\_\{2\},\\dots,\\lVert\\bm\{\\theta\}\_\{c,\\mathrm\{tr\}\}\\rVert\_\{2\}\)andΘ~tr=Θtr​S−1\\widetilde\{\\Theta\}\_\{\\mathrm\{tr\}\}=\\Theta\_\{\\mathrm\{tr\}\}S^\{\-1\}, the inner column\-normalised ridge fit is

𝝃^tr=S−1​\(Θ~tr⊤​Θ~tr\+λ​I\)−1​Θ~tr⊤​𝒃tr\.\\widehat\{\\bm\{\\xi\}\}\_\{\\mathrm\{tr\}\}=S^\{\-1\}\\bigl\(\\widetilde\{\\Theta\}\_\{\\mathrm\{tr\}\}^\{\\\!\\top\}\\widetilde\{\\Theta\}\_\{\\mathrm\{tr\}\}\+\\lambda I\\bigr\)^\{\-1\}\\widetilde\{\\Theta\}\_\{\\mathrm\{tr\}\}^\{\\\!\\top\}\\bm\{b\}\_\{\\mathrm\{tr\}\}\.\(20\)Appendix[9\.1](https://arxiv.org/html/2608.12879#S9.SS1)gives the row split, normalisations, duplicate penalty, and numerical constants\. Differential evolution optimises the continuous orders within each temporal branch[20](https://arxiv.org/html/2608.12879#bib.bib14); in the integer branch,α=1\\alpha=1is fixed\. Powers are enumerated only as nondecreasing patterns𝔓c\\mathfrak\{P\}\_\{c\}, eliminating permutation duplicates\. This bi\-level design confines the global search to a small order space while solving the linear coefficients at every objective evaluation\.

#### 4\.3The parsimony\-promoting elbow

Solving \([19](https://arxiv.org/html/2608.12879#S4.E19)\) for increasingccproduces a sequence of best modelsℳ1⋆,ℳ2⋆,…\\mathcal\{M\}\_\{1\}^\{\\star\},\\mathcal\{M\}\_\{2\}^\{\\star\},\\dotsand a Pareto front of validation error against support size\. Model quality is scored on held\-out evaluation rows,

ℰval​\(ℳ\)=nval−1​∥𝒃val−Θval​𝝃^∥22Var⁡\(𝒃val\)\+ε,\\mathcal\{E\}\_\{\\mathrm\{val\}\}\(\\mathcal\{M\}\)=\\frac\{n\_\{\\mathrm\{val\}\}^\{\-1\}\\lVert\\bm\{b\}\_\{\\mathrm\{val\}\}\-\\Theta\_\{\\mathrm\{val\}\}\\widehat\{\\bm\{\\xi\}\}\\rVert\_\{2\}^\{2\}\}\{\\operatorname\{Var\}\(\\bm\{b\}\_\{\\mathrm\{val\}\}\)\+\\varepsilon\},\(21\)whereVar⁡\(𝒃val\)=nval−1​∑i=1nval\(bi−b¯val\)2\\operatorname\{Var\}\(\\bm\{b\}\_\{\\mathrm\{val\}\}\)=n\_\{\\mathrm\{val\}\}^\{\-1\}\\sum\_\{i=1\}^\{n\_\{\\mathrm\{val\}\}\}\(b\_\{i\}\-\\overline\{b\}\_\{\\mathrm\{val\}\}\)^\{2\}is the empirical validation\-target variance used by the implementation, andε=10−14\\varepsilon=10^\{\-14\}guards against a zero or numerically negligible denominator\. Variance normalisation makes the score dimensionless and comparable across temporal targets whose scale changes withα\\alpha\. Training and validation use disjoint weak rows, although overlapping test\-function supports make them correlated measurements of the same noisy field\. ThusKKcounts weak equations and is distinct from an effective number of independent observations\. We useℰval\\mathcal\{E\}\_\{\\mathrm\{val\}\}for internal model selection and multi\-seed experiments for across\-realisation assessment\.

Letℰc⋆\\mathcal\{E\}\_\{c\}^\{\\star\}be the best validation error at support sizecc\. The sweep stops when the relative improvement falls belowδ\\deltaafter the earliest admissible stopping sizecmin=2c\_\{\\min\}=2, or when the selected elbow remains unchanged after an additional size\. To define the elbow, the implementation uses the benefit coordinateyc=−log10⁡\(ℰc⋆\+ε\)y\_\{c\}=\-\\log\_\{10\}\(\\mathcal\{E\}\_\{c\}^\{\\star\}\+\\varepsilon\)together with support sizecc\. After min–max normalising both coordinates, the selected interior point maximises its signed vertical excess above the chord joining the first and last points\. This chord\-based criterion is related to normalised difference\-curve knee detection[15](https://arxiv.org/html/2608.12879#bib.bib20)\. Thus “above” refers to the normalised benefit–complexity coordinates; when the same front is drawn as raw validation error versus support size, the corresponding knee appears below the endpoint chord\. If no interior point has positive signed excess, the smallest model is retained\. With only two sizes, the larger model is selected only when it improveslog10⁡ℰc⋆\\log\_\{10\}\\mathcal\{E\}\_\{c\}^\{\\star\}by at least0\.150\.15\. The two\-point margin is heuristic; Appendix[11\.1](https://arxiv.org/html/2608.12879#S11.SS1)examines sensitivity over0\.100\.10–0\.200\.20\.

Algorithm 1Weak FPDE discovery via Pareto\-based subset selection1:Data

uuon a grid; test functions

\{ϕk\}\\\{\\phi\_\{k\}\\\}; powers

𝒜=\{0,…,P\}\\mathcal\{A\}=\\\{0,\\dots,P\\\}, which define

𝔓c\\mathfrak\{P\}\_\{c\}; max size

cmaxc\_\{\\max\}; earliest stopping size

cmin=2c\_\{\\min\}=2; plateau tolerance

δ\\delta
2:Build the weak target

𝒃⁡\(mα,α\)\\bm\{b\}\(m\_\{\\alpha\},\\alpha\)and weak columns

Θk,j\\Theta\_\{k,j\}via the adjoints \([9](https://arxiv.org/html/2608.12879#S3.E9)\)–\([11](https://arxiv.org/html/2608.12879#S3.E11)\)

3:Split rows into training/validation

4:for

c=1,2,…,cmaxc=1,2,\\dots,c\_\{\\max\}do

5:foreach

𝒑∈𝔓c\\bm\{p\}\\in\\mathfrak\{P\}\_\{c\}generated from

𝒜\\mathcal\{A\}do

6:foreach admitted temporal mode

mαm\_\{\\alpha\}do

7:Optimise

α\\alphaonly within

ℐmα\\mathcal\{I\}\_\{m\_\{\\alpha\}\}\(fixed at

11for

mα=intm\_\{\\alpha\}=\\mathrm\{int\}\) and optimise

𝜷\\bm\{\\beta\}, refitting

𝝃^\\widehat\{\\bm\{\\xi\}\}by \([20](https://arxiv.org/html/2608.12879#S4.E20)\)

8:endfor

9:Retain the lowest\-objective temporal mode for this power pattern

10:endfor

11:

ℳc⋆←\\mathcal\{M\}\_\{c\}^\{\\star\}\\leftarrowbest model of size

ccby the objective

JcJ\_\{c\}in \([19](https://arxiv.org/html/2608.12879#S4.E19)\); let

ℰc⋆=ℰval​\(ℳc⋆\)\\mathcal\{E\}^\{\\star\}\_\{c\}=\\mathcal\{E\}\_\{\\mathrm\{val\}\}\(\\mathcal\{M\}\_\{c\}^\{\\star\}\)
12:Recompute the current signed\-chord elbow

c^elbow\\widehat\{c\}\_\{\\mathrm\{elbow\}\}from

\{ℳj⋆\}j=1c\\\{\\mathcal\{M\}\_\{j\}^\{\\star\}\\\}\_\{j=1\}^\{c\}
13:if

c≥cminc\\geq c\_\{\\min\}and the validation\-improvement plateau condition holdsthen

14:break

15:elseif

c^elbow<c\\widehat\{c\}\_\{\\mathrm\{elbow\}\}<cand the elbow is unchanged for the declared patiencethen

16:break⊳\\trianglerightselection stability

17:endif

18:endfor

19:Select the elbow model

ℳ⋆\\mathcal\{M\}^\{\\star\}of the Pareto front

\{ℳc⋆\}\\\{\\mathcal\{M\}\_\{c\}^\{\\star\}\\\}
20:Prune numerically inactive terms \(Section[4\.4](https://arxiv.org/html/2608.12879#S4.SS4)\) and refit at exact orders \(Section[4\.5](https://arxiv.org/html/2608.12879#S4.SS5)\)

21:return

ℳ⋆\\mathcal\{M\}^\{\\star\}⊳\\trianglerightrepresenting the discovered fractional equation

#### 4\.4Inactive\-term pruning

A selected model can occasionally contain a term with negligible fitted effect\. We remove such terms using a scale\-aware, non\-oracle rule\. ForΘ​𝝃^=∑jξ^j​𝜽j\\Theta\\widehat\{\\bm\{\\xi\}\}=\\sum\_\{j\}\\widehat\{\\xi\}\_\{j\}\\bm\{\\theta\}\_\{j\}, definerj=∥ξ^j​𝜽j∥2/\(∥Θ​𝝃^∥2\+ε\)r\_\{j\}=\\lVert\\widehat\{\\xi\}\_\{j\}\\bm\{\\theta\}\_\{j\}\\rVert\_\{2\}/\(\\lVert\\Theta\\widehat\{\\bm\{\\xi\}\}\\rVert\_\{2\}\+\\varepsilon\)\. Termjjis removed whenrj≤τcontribr\_\{j\}\\leq\\tau\_\{\\mathrm\{contrib\}\}or\|ξ^j\|≤τabs\|\\widehat\{\\xi\}\_\{j\}\|\\leq\\tau\_\{\\mathrm\{abs\}\}\. The ratio is evaluated over all finite weak rows\. Unlike raw coefficient magnitude, it accounts for the different scales of fractional\-library columns\. Because partially cancelling terms can makerj\>1r\_\{j\}\>1, only small values are interpreted\.

#### 4\.5Exact\-order refitting

During the global search, weak columns are interpolated from features precomputed on an order grid\. Here “exact\-order” means that the operators are subsequently evaluated directly at continuous order values instead of by interpolation\. After elbow selection and pruning, local gradient\-free steps refine the selected orders within the temporal branch chosen during the global search\. If a fractional branch is selected,α\\alphaand the non\-identity spatial orders are polished within their branch bounds and local trust regions\. If the exact\-integer branch is selected,α\\alpharemains fixed at11and only the non\-identity spatial orders are polished\. The temporal\-branch comparison is completed during model selection, so this conditional refit operates only within the selected branch\. Letℳ⋆\\mathcal\{M\}^\{\\star\}denote the selected model and\(α^ℳ⋆,𝜷^ℳ⋆\)\(\\widehat\{\\alpha\}^\{\\,\\mathcal\{M\}^\{\\star\}\},\\widehat\{\\bm\{\\beta\}\}^\{\\,\\mathcal\{M\}^\{\\star\}\}\)its polished orders\. The final coefficients are then

𝝃^ℳ⋆=arg​min𝝃∥𝐛ℳ⋆−Θℳ⋆𝝃∥22\+λ∥Sℳ⋆𝝃∥22,\\widehat\{\\bm\{\\xi\}\}^\{\\,\\mathcal\{M\}^\{\\star\}\}=\\argmin\_\{\\bm\{\\xi\}\}\\;\\bigl\\lVert\\bm\{b\}^\{\\mathcal\{M\}^\{\\star\}\}\-\\Theta^\{\\mathcal\{M\}^\{\\star\}\}\\bm\{\\xi\}\\bigr\\rVert\_\{2\}^\{2\}\+\\lambda\\bigl\\lVert S^\{\\mathcal\{M\}^\{\\star\}\}\\bm\{\\xi\}\\bigr\\rVert\_\{2\}^\{2\},\(22\)where the target and design are evaluated directly at the polished orders andSℳ⋆S^\{\\mathcal\{M\}^\{\\star\}\}column\-normalises the design\. The same ridge parameter is used during search and refitting\. Model selection is completed before this update; the training score, validation score, objective, and heuristic information criteria therefore keep their selection\-stage meanings\. The final all\-row residual is stored separately asℰfit\\mathcal\{E\}\_\{\\mathrm\{fit\}\}\. Direct evaluation removes order\-grid interpolation error from this conditional refit\. The selected active\-term set, temporal mode, and optimisation basin can still depend on the order grid, search trajectory, order bounds, test functions, and data resolution\.

#### 4\.6Computational cost

Two choices control the cost\. First, the bi\-level formulation restricts global optimisation to the order variables—c\+1c\+1dimensions for a fractional temporal branch andccfor the exact integer branch—while the coefficients are solved directly\. Second, early stopping usually evaluates support sizes only up to the selected model plus one\.

LetN=nt​nxN=n\_\{t\}n\_\{x\}be the number of field samples,K=Kt​KxK=K\_\{t\}K\_\{x\}the number of weak rows,GαG\_\{\\alpha\}andGβG\_\{\\beta\}the numbers of temporal and spatial order nodes, andP\+1P\+1the number of powers\. Weak\-Pareto precomputes the target and candidate features at these order nodes\. At support sizecc, the number of non\-redundant power patterns isNp​\(c\)=\(P\+cc\)N\_\{p\}\(c\)=\\binom\{P\+c\}\{c\}\. A separable projectionT​U​X⊤TUX^\{\\\!\\top\}costsO⁡\(Kt​N\+K​nx\)O\(K\_\{t\}N\+Kn\_\{x\}\), and a periodic spectral adjoint at one order costsO⁡\(Kx​nx​log⁡nx\)O\(K\_\{x\}n\_\{x\}\\log n\_\{x\}\)\. Precomputation is therefore linear inGαG\_\{\\alpha\}andGβG\_\{\\beta\}, with memoryO⁡\(K⁡\[Gα\+\(P\+1\)​Gβ\]\+N\)O\\\!\\left\(K\[G\_\{\\alpha\}\+\(P\+1\)G\_\{\\beta\}\]\+N\\right\)\.

After precomputation, one objective evaluation forms aK×cK\\times cdesign and solves a ridge system inO⁡\(K​c2\+c3\)O\(Kc^\{2\}\+c^\{3\}\)time andO⁡\(K​c\)O\(Kc\)memory\. IfNDE​\(c\)N\_\{\\mathrm\{DE\}\}\(c\)is the number of differential\-evolution evaluations per power pattern and the sweep stops atcstopc\_\{\\mathrm\{stop\}\}, the total search cost is

O⁡\(∑c=1cstopNp​\(c\)​NDE​\(c\)​\[K​c2\+c3\]\)O\\\!\\left\(\\sum\_\{c=1\}^\{c\_\{\\mathrm\{stop\}\}\}N\_\{p\}\(c\)N\_\{\\mathrm\{DE\}\}\(c\)\\,\[Kc^\{2\}\+c^\{3\}\]\\right\)after precomputation\. In the reported experimentsc≤4c\\leq 4; hence the cubic term is negligible and each objective evaluation is effectively linear in the number of weak rows\. Fast Fourier transforms accelerate the dominant periodic\-operator calculations\.

### 5Experiments and results

The experiments address three primary questions\. First, can Weak\-Pareto recover parsimonious FPDEs across linear and nonlinear benchmarks? Second, are weak measurements more robust than pointwise fractional features under matched selection? Third, is continuous\-order subset search more reliable than a dense fixed dictionary? Supporting studies examine the superunit temporal branch, a two\-dimensional anisotropic example, a contemporary neural fractional\-discovery framework, computational cost, and applicability to irregular experimental data\.

#### 5\.1Empirical benchmarks and evaluation metrics

We use four periodic main benchmarks \(Table[2](https://arxiv.org/html/2608.12879#S5.T2)\): FADE, two Riesz reaction–diffusion \(RD\) equations with integer or fractional time, and a nonlinear fractional Burgers equation\. The Burgers case is the principal nonlinear benchmark because its support contains the genuine quadratic transport term−u∂xu\-u\\,\\partial\_\{x\}u; it also tests discrimination between fractional diffusion of order1\.71\.7and the nearby integer second derivative\. Section[5\.3](https://arxiv.org/html/2608.12879#S5.SS3)reports a separate fixed\-support superunit diagnostic, and Appendix[11](https://arxiv.org/html/2608.12879#S11)adds an integer\-only equation and a case with two fractional spatial derivatives\. Unless stated otherwise, noise is multiplicative,u↦u⁡\(1\+ρ​ζ\)u\\mapsto u\(1\+\\rho\\zeta\)withζ∼𝒰⁡\[−1,1\]\\zeta\\sim\\mathcal\{U\}\[\-1,1\], and tables report100​ρ100\\rhopercent\. Recovery counts are over five seeds, order and coefficient errors are conditioned on support/power recovery, andℰfit\\mathcal\{E\}\_\{\\mathrm\{fit\}\}is summarised over all seeds\. The10%10\\%Riesz cases are deliberately severe identifiability tests at the present resolution\. The time–space reaction–diffusion generator and evaluator share the Caputo L1 discretisation; therefore, the clean experiment is a discretisation\-consistency test rather than an independent forward\-solver validation\.

The main experiments use one common, non\-oracle configuration: automatic support\-size stopping,cmax=4c\_\{\\max\}=4, candidate powersp∈\{0,1,2\}p\\in\\\{0,1,2\\\}, and elbow selection\. Appendix[11](https://arxiv.org/html/2608.12879#S11)examines sensitivity to these choices\.

##### Evaluation protocol\.

We separate structural recovery from parameter accuracy\.*Support/power recovery*requires the correct number of terms and matching integer powers after pruning\.*Operator\-structure recovery*additionally requires the correct temporal branch, spatial operator modes, and an absolute error no greater thanτq=0\.15\\tau\_\{q\}=0\.15for every positive derivative order; identity terms must be identified exactly\. A fractional time order near one is not the exact operator∂t\\partial\_\{t\}, and a low\-order Riesz term is not the reaction identity\.

The tolerance is a predefined operational criterion\. For the main\-text synthetic benchmarks, the positive true orders range from approximately0\.80\.8to2\.02\.0; henceτq=0\.15\\tau\_\{q\}=0\.15corresponds to relative deviations of7\.5%7\.5\\%–18\.75%18\.75\\%\. An absolute criterion is appropriate because fractional order is dimensionless and temporal and spatial orders are measured on the same scale\. Rescoring the 60 noisy runs in the matched weak\-versus\-strong comparisons gives 45, 47, and 47 recoveries for Weak\-Pareto atτq=0\.125\\tau\_\{q\}=0\.125,0\.150\.15, and0\.1750\.175, respectively; the strong\-form method gives 0/60 at all three values\. Thus, the main comparison is insensitive to moderate changes around the reported tolerance\. We report continuous order errors alongside the binary counts because the cut\-off is an interpretation rule, not a measure of uncertainty\.

Parameter accuracy is measured byeα=\|α^−α⋆\|e\_\{\\alpha\}=\|\\widehat\{\\alpha\}\-\\alpha^\{\\star\}\|\. To compare right\-hand\-side terms, we define a one\-to\-one matching mapπ\\pifrom true\-term indices to selected\-term indices\. Processing the true terms in their stored order,π⁡\(j\)\\pi\(j\)assigns termjjto the nearest unmatched selected order with the same power\. With𝝃^π=\(ξ^π⁡\(1\),…,ξ^π⁡\(c\)\)\\widehat\{\\bm\{\\xi\}\}\_\{\\pi\}=\(\\widehat\{\\xi\}\_\{\\pi\(1\)\},\\dots,\\widehat\{\\xi\}\_\{\\pi\(c\)\}\)andεξ=10−12\\varepsilon\_\{\\xi\}=10^\{\-12\}, we use

eβmax=maxj⁡\|β^π⁡\(j\)−βj⋆\|,eξmax=maxj⁡\|ξ^π⁡\(j\)−ξj⋆\|\|ξj⋆\|\+εξ,eξ,2=∥𝝃^π−𝝃⋆∥2∥𝝃⋆∥2\+εξ\.e\_\{\\beta\}^\{\\max\}=\\max\_\{j\}\|\\widehat\{\\beta\}\_\{\\pi\(j\)\}\-\\beta\_\{j\}^\{\\star\}\|,\\hskip 20\.00003pte\_\{\\xi\}^\{\\max\}=\\max\_\{j\}\\frac\{\|\\widehat\{\\xi\}\_\{\\pi\(j\)\}\-\\xi\_\{j\}^\{\\star\}\|\}\{\|\\xi\_\{j\}^\{\\star\}\|\+\\varepsilon\_\{\\xi\}\},\\hskip 20\.00003pte\_\{\\xi,2\}=\\frac\{\\lVert\\widehat\{\\bm\{\\xi\}\}\_\{\\pi\}\-\\bm\{\\xi\}^\{\\star\}\\rVert\_\{2\}\}\{\\lVert\\bm\{\\xi\}^\{\\star\}\\rVert\_\{2\}\+\\varepsilon\_\{\\xi\}\}\.These errors are averaged only over runs with correct support and powers;eξ,2e\_\{\\xi,2\}is especially useful when a true coefficient is small\. Model selection uses the held\-out, variance\-normalised scoreℰval\\mathcal\{E\}\_\{\\mathrm\{val\}\}\. After selection, the model is refit on all weak rows and reported withℰfit=∥𝒃−Θ​𝝃^∥2/\(∥𝒃∥2\+ε\)\\mathcal\{E\}\_\{\\mathrm\{fit\}\}=\\lVert\\bm\{b\}\-\\Theta\\widehat\{\\bm\{\\xi\}\}\\rVert\_\{2\}/\(\\lVert\\bm\{b\}\\rVert\_\{2\}\+\\varepsilon\)\. Because the weak and strong frameworks use different regression rows,ℰfit\\mathcal\{E\}\_\{\\mathrm\{fit\}\}is a within\-framework diagnostic; cross\-method conclusions rely on recovery rates and parameter errors\.

Table 2:Main benchmark equations in the encoding of Eq\. \([4](https://arxiv.org/html/2608.12879#S2.E4)\);ℛβ\\mathcal\{R\}\_\{\\beta\}is the Riesz operator and\(0,0\)\(0,0\)denotes the identity termuu\. The spatial operator is directional for FADE and fractional Burgers and Riesz for the two reaction–diffusion benchmarks\. True terms are listed directly as triples\(p,β,ξ\)\(p,\\beta,\\xi\); horizons\(T,Lx\)\(T,L\_\{x\}\)are approximately\(15,30\)\(15,30\),\(1,20\)\(1,20\),\(3,20\)\(3,20\), and\(12,30\)\(12,30\), respectively\.

#### 5\.2Recovery accuracy of Weak\-Pareto

Table[3](https://arxiv.org/html/2608.12879#S5.T3)reports Weak\-Pareto at10%10\\%multiplicative noise, used as a common stress level across the four benchmarks\. FADE and fractional Burgers are recovered in all five seeds, including the operator orders, with small spatial\-order and coefficient\-vector errors\. The Riesz reaction–diffusion cases are more difficult\. Support and powers are recovered in3/53/5space\-fractional runs and all five time–space runs, but no run satisfies the full operator criterion at10%10\\%noise\. Order identification is the more persistent difficulty, although support recovery also degrades in the space\-fractional case at10%10\\%noise; the worst relative coefficient error is amplified further by the small reaction coefficients\. Table[4](https://arxiv.org/html/2608.12879#S5.T4)resolves these effects across lower noise levels\.

Table 3:Weak\-Pareto on the four main benchmarks at10%10\\%multiplicative noise\. Recovery definitions and reporting conventions follow the evaluation protocol in Section[5\.1](https://arxiv.org/html/2608.12879#S5.SS1)\. The relativeeξmaxe\_\{\\xi\}^\{\\max\}is inflated on the reaction–diffusion benchmarks by their small reaction coefficient, for whicheξ,2e\_\{\\xi,2\}is more informative\. Temporal\-bound concentration is detailed in Table[4](https://arxiv.org/html/2608.12879#S5.T4)\.Table[4](https://arxiv.org/html/2608.12879#S5.T4)now includes the clean reference for both Riesz reaction–diffusion benchmarks\. At0%0\\%noise, both recover the complete operator in all five seeds, confirming that the subsequent failures are noise\-induced\. The clean time–spaceeαe\_\{\\alpha\}standard deviation is1\.4×10−61\.4\\times 10^\{\-6\}before rounding, although it appears as0\.00000\.0000at the precision used in Table[4](https://arxiv.org/html/2608.12879#S5.T4)\. Lowering the positive noise level improves the Riesz\-order estimates, but complete operator recovery remains difficult because it also requires the correct temporal mode and the reaction identity\. The time–space case reaches2/52/5operator recoveries at2%2\\%noise; no positive\-noise row for the space\-fractional case does\. Boundary attainment is clearest for the space\-fractional case: all five support/power\-recovered runs at5%5\\%noise and all three runs entering the conditioned10%10\\%summary haveeα=0\.20±0\.00e\_\{\\alpha\}=0\.20\\pm 0\.00\. Since the true order isα=1\.00\\alpha=1\.00andαmin=0\.80\\alpha\_\{\\min\}=0\.80, this is the truncation distance; zero dispersion therefore indicates constraint saturation, not precision\. In the time–space case,eαe\_\{\\alpha\}rises from0\.07±0\.030\.07\\pm 0\.03at2%2\\%noise to0\.17±0\.010\.17\\pm 0\.01at10%10\\%, indicating increasing concentration near the lower boundary\. Order errors are conditioned on support/power recovery\. Appendix[8](https://arxiv.org/html/2608.12879#S8)identifies the noisy temporal target as the dominant source of the shift in these profiles\. The positive\-noise rows should therefore be read as support recovery with increasingly accurate spatial order at lower noise\.

Table 4:Order identifiability of the two Riesz reaction–diffusion benchmarks versus noise\. Reporting conventions follow the evaluation protocol in Section[5\.1](https://arxiv.org/html/2608.12879#S5.SS1)\.##### Support\-conditioned spatial\-order diagnostic\.

We next test whether the local sensitivity in Eq\. \([16](https://arxiv.org/html/2608.12879#S3.E16)\) is consistent with the spatial\-order errors above\. For each main benchmark,SβS\_\{\\beta\}is evaluated on the clean field and true support using the exact weak features; the temporal branch and all other orders are fixed at their true values\. In a separate fixed\-support profile under multiplicative noise, only the principal noninteger linear spatial order is varied, while all coefficients are refitted on the same training rows with the ridge parameter and variance\-normalised validation score used by Weak\-Pareto\. This removes support\-selection and temporal\-order errors from the diagnostic\. Table[5](https://arxiv.org/html/2608.12879#S5.T5)reports the results\.

Table 5:Support\-conditioned local sensitivitySβS\_\{\\beta\}and fixed\-support absolute erroreβ,fixe\_\{\\beta,\\mathrm\{fix\}\}of the principal noninteger linear spatial order \(five seeds\)\. The complete0%0\\%,2%2\\%,5%5\\%, and10%10\\%profiles and per\-seed estimates are provided in Online Resource 1\.The sensitivity ordering is the reverse of the noisy fixed\-support error ordering: fractional Burgers has the largestSβS\_\{\\beta\}and smallest10%10\\%error, followed by FADE, whereas the two Riesz benchmarks have much smaller sensitivities and substantially larger errors\. The same qualitative ordering appears in the complete\-discoveryeβmaxe\_\{\\beta\}^\{\\max\}values in Tables[3](https://arxiv.org/html/2608.12879#S5.T3)and[4](https://arxiv.org/html/2608.12879#S5.T4)\. All four fixed\-support profiles remain accurate without noise\. The four\-benchmark ordering therefore provides an explanatory consistency check:SβS\_\{\\beta\}characterises local, support\-conditioned susceptibility to spatial\-order perturbations\. Caputo endpoint sensitivity and dense\-dictionary collinearity arise from separate mechanisms analysed elsewhere\.

#### 5\.3Superunit temporal\-order diagnostic

The main benchmark suite contains no true temporal order in\(1,2\)\(1,2\)\. We therefore test the superunit branch on semi\-analytic periodic data satisfying

Dt1\.650C​u=0\.12​Dx2​u,∂tu⁡\(0,x\)=0\.\{\}\_\{0\}^\{C\}D\_\{t\}^\{1\.65\}u=0\.12\\,D\_\{x\}^\{2\}u,\\hskip 20\.00003pt\\partial\_\{t\}u\(0,x\)=0\.\(23\)For a spatial Fourier mode with wavenumberκ\\kappa, the mode amplitude evolves asE1\.65,1​\(−0\.12​κ2​t1\.65\)E\_\{1\.65,1\}\(\-0\.12\\kappa^\{2\}t^\{1\.65\}\), whereEa,b​\(z\)=∑m=0∞zm/Γ⁡\(a​m\+b\)E\_\{a,b\}\(z\)=\\sum\_\{m=0\}^\{\\infty\}z^\{m\}/\\Gamma\(am\+b\)is the two\-parameter Mittag–Leffler function\. This semi\-analytic evolution is independent of the L1 temporal discretisation used by the discovery evaluator\. To isolate temporal\-branch and order recovery, the support is fixed to one linear term \(c=1c=1,p1=0p\_\{1\}=0\), while both methods search the temporal branch,α\\alpha, andβ\\betaand fit the coefficient\. This is therefore a branch/order diagnostic\. We use the same noisy realisation for the weak and strong frameworks at each of five seeds\. Operator recovery requires the superunit branch together with\|α^−1\.65\|≤τq\|\\widehat\{\\alpha\}\-1\.65\|\\leq\\tau\_\{q\}and\|β^−2\|≤τq\|\\widehat\{\\beta\}\-2\|\\leq\\tau\_\{q\}\. Table[6](https://arxiv.org/html/2608.12879#S5.T6)reports the resulting branch and operator recoveries\.

Table 6:Fixed\-support superunit diagnostic for Eq\. \([23](https://arxiv.org/html/2608.12879#S5.E23)\) under multiplicative\-uniform noise\. Complete per\-seed estimates and mean±\\pmsample\-standard\-deviation errors are provided in Online Resource 1\.Both methods recover the clean superunit operator in all five runs; for Weak\-Pareto, the clean temporal\- and spatial\-order errors are below10−310^\{\-3\}and the mean relative coefficient error is0\.0280\.028\. At0\.5%0\.5\\%noise, Weak\-Pareto retains 5/5 branch and operator recovery, witheα=0\.074±0\.038e\_\{\\alpha\}=0\.074\\pm 0\.038,eβ=0\.037±0\.028e\_\{\\beta\}=0\.037\\pm 0\.028, andeξ=0\.073±0\.064e\_\{\\xi\}=0\.073\\pm 0\.064, whereas the strong\-form comparator selects the wrong temporal branch in every run\. At1%1\\%, Weak\-Pareto still selects the superunit branch in all five runs, but four estimates reach the upper search boundα=1\.85\\alpha=1\.85, leaving only 1/5 complete operator recoveries andeα=0\.172±0\.062e\_\{\\alpha\}=0\.172\\pm 0\.062\. The diagnostic therefore verifies that the branch\-aware weak search extends toα\>1\\alpha\>1under clean and mild noise\. The upper\-bound saturation at1%1\\%is consistent with the endpoint sensitivity analysed in Remark[2](https://arxiv.org/html/2608.12879#Thmremark2)\.

#### 5\.4Robustness relative to a strong\-form library

Fig\.[2](https://arxiv.org/html/2608.12879#S5.F2)and Table[7](https://arxiv.org/html/2608.12879#S5.T7)compare Weak\-Pareto with a strong\-form fractional library under the same best\-subset framework on FADE and fractional Burgers\. The methods are comparable without noise\. Under multiplicative noise, however, Weak\-Pareto recovers the correct support in all five seeds at every tested level up to20%20\\%on both benchmarks\. The strong\-form framework does not recover the correct Burgers support in any noisy run and in all but the5%5\\%FADE condition, where it recovers4/54/5supports but has order\-one coefficient error \(eξmax≈1\.1e\_\{\\xi\}^\{\\max\}\\approx 1\.1\)\. Weak\-Pareto’s Burgers coefficient error remains below0\.0220\.022throughout; on FADE it remains small through10%10\\%noise and rises only at20%20\\%, when the smaller diffusion coefficient becomes difficult to estimate\.

The complete\-framework comparison includes Weak\-Pareto’s exact\-order refinement\. The library\-only contrast in Section[5\.10](https://arxiv.org/html/2608.12879#S5.SS10)instead compares Strong\-Pareto with Weak\-Pareto without polishing under the same selector, isolating pointwise versus weak candidate measurements\. Operator recovery then changes from0/50/5to5/55/5at10%10\\%FADE noise, identifying the weak library as the primary source of robustness\. Appendix[13](https://arxiv.org/html/2608.12879#S13)reaches the same conclusion under additive Gaussian noise: Weak\-Pareto recovers the correct support in all five runs and four of five complete operators on both benchmarks, while the strong\-form framework achieves no correct\-support recovery\.

Appendix[17](https://arxiv.org/html/2608.12879#S17)provides a complementary comparison with the adapted neural fractional\-discovery framework of Yu et al\.[23](https://arxiv.org/html/2608.12879#bib.bib19)on the advection–diffusion benchmark\. We treat it as a method\-level comparison because the two methods use different fractional\-operator realisations; the matched weak–strong experiment provides the operator\-controlled ablation\.

Figure 2:Relative coefficient erroreξmaxe\_\{\\xi\}^\{\\max\}versus noise on \(a\) FADE and \(b\) fractional Burgers for the weak and strong\-form libraries under the same selector; lower is better\. The solid curve shows Weak\-Pareto\. Dotted horizontal lines show the strong\-form result at0%0\\%noise in both panels and its additional recovered case on FADE at5%5\\%noiseTable 7:Weak\-Pareto versus the strong\-form pointwise library under the same best\-subset selector on FADE\. The fit residual is a within\-framework diagnostic\.The matched library\-only ablation discussed in Section[5\.10](https://arxiv.org/html/2608.12879#S5.SS10)provides the direct evidence for the weak formulation: with the selector held fixed, replacing pointwise features by weak measurements changes FADE operator recovery at10%10\\%noise from0/50/5to5/55/5\. The additive\-Gaussian experiment in Appendix[13](https://arxiv.org/html/2608.12879#S13)confirms that this advantage is not tied to the multiplicative\-uniform noise model\.

#### 5\.5Model selection via Pareto\-based subset selection

Fig\.[3](https://arxiv.org/html/2608.12879#S5.F3)and Table[8](https://arxiv.org/html/2608.12879#S5.T8)illustrate the support\-size search on FADE at5%5\\%noise\. Validation error drops by more than an order of magnitude from one to the true two\-term support and improves only marginally at three terms\. The elbow therefore selectsc=2c=2, and selection stability stops the search after evaluatingc=3c=3, before the permitted maximumcmax=4c\_\{\\max\}=4\. The procedure thus expresses parsimony through an explicit, data\-driven stopping rule\.

Figure 3:Validation errorℰval\\mathcal\{E\}\_\{\\mathrm\{val\}\}versus support size on FADE \(5%5\\%noise\); the marked elbow \(c=2c=2\) is the selected model\. The search halts automatically afterc=3c=3; larger models up tocmax=4c\_\{\\max\}=4are therefore not exploredTable 8:Support\-size progress on FADE \(5%5\\%noise\): training and validation error of the best model at each support size\. The signed elbow selectsc=2c=2\.
#### 5\.6A nonlinear fractional benchmark

The fractional Burgers equation∂tu=−u∂xu\+0\.25Dx1\.7u\\partial\_\{t\}u=\-u\\,\\partial\_\{x\}u\+0\.25D\_\{x\}^\{1\.7\}uis the principal test of nonlinear discovery because its support contains the genuine quadratic transport term−u∂xu\-u\\,\\partial\_\{x\}u\. It also tests whether the framework distinguishes the fractional diffusion termDx1\.7​uD\_\{x\}^\{1\.7\}ufrom plausible integer\-order and nonlinear alternatives\. We compare the true two\-term structure with competing two\-term structures assembled fromuxu\_\{x\},u2​uxu^\{2\}u\_\{x\},ux​xu\_\{xx\},u​D1\.7​uuD^\{1\.7\}u, andD0\.5​uD^\{0\.5\}u\. For each structure, coefficients are refit at the exact candidate orders and scored by the relative weak residual\. The residual margin is the smallest residual among the structures that do not match the ground truth, divided by the true\-structure residual\. Weak\-Pareto selects the true pair at all six tested noise levels \(0%0\\%,5%5\\%,10%10\\%,15%15\\%,20%20\\%, and25%25\\%\)\. The margin decreases monotonically from218×218\\timeswithout noise to3\.2×3\.2\\timesat25%25\\%noise but remains above one, while the coefficient error stays below0\.020\.02\. The closest competing structure replaces fractional diffusion byux​xu\_\{xx\}at every level, showing that the weak library continues to distinguish both the noninteger diffusion order and the nonlinear transport term under substantial noise\. Fig\.[4](https://arxiv.org/html/2608.12879#S5.F4)shows the complete six\-point residual curve, whereas Table[9](https://arxiv.org/html/2608.12879#S5.T9)reports the representative0%0\\%,10%10\\%, and25%25\\%rows together with coefficient errors\.

Table 9:Representative noise levels for nonlinear fractional Burgers\. The residual margin is the smallest residual among the structures that do not match the ground truth, divided by the true\-structure residual;eξmaxe\_\{\\xi\}^\{\\max\}is the relative coefficient error of the true structure\.![Refer to caption](https://arxiv.org/html/2608.12879v1/Fig4.png)Figure 4:\(a\) Space–time solution field of the nonlinear fractional Burgers benchmark\. \(b\) Weak residual of the true two\-term structure and the closest competing structure at0%0\\%,5%5\\%,10%10\\%,15%15\\%,20%20\\%, and25%25\\%noise on a logarithmic scale; the true structure remains separated throughout
#### 5\.7Experimental frozen\-soil creep

We next test the weak\-form principle on naturally noisy, irregularly sampled clay and silt creep measurements from Yu et al\.[23](https://arxiv.org/html/2608.12879#bib.bib19)\. After averaging one repeated silt timestamp, the records contain 56 and 43 strain observations under constant loads of1\.111\.11and1\.141\.14MPa, respectively\. We use the fractional Kelvin model

DαtC​ϵ​\(t\)=σloadηK−EηK​ϵ​\(t\),0<α<1,\{\}^\{C\}D\_\{t\}^\{\\alpha\}\\epsilon\(t\)=\\frac\{\\sigma\_\{\\mathrm\{load\}\}\}\{\\eta\_\{\\mathrm\{K\}\}\}\-\\frac\{E\}\{\\eta\_\{\\mathrm\{K\}\}\}\\epsilon\(t\),\\hskip 20\.00003pt0<\\alpha<1,\(24\)whereϵ\\epsilonis strain,EEis the elastic modulus,ηK\\eta\_\{\\mathrm\{K\}\}is the Kelvin viscosity parameter, andσload\\sigma\_\{\\mathrm\{load\}\}is the applied constant stress\. Applying the fractional integral gives the smoothing representation

ϵ⁡\(t\)−ϵ⁡\(0\)=a0​tαΓ⁡\(α\+1\)\+a1​Itα​ϵ​\(t\),a0=σloadηK,a1=−EηK\.\\epsilon\(t\)\-\\epsilon\(0\)=a\_\{0\}\\frac\{t^\{\\alpha\}\}\{\\Gamma\(\\alpha\+1\)\}\+a\_\{1\}I\_\{t\}^\{\\alpha\}\\epsilon\(t\),\\hskip 20\.00003pta\_\{0\}=\\frac\{\\sigma\_\{\\mathrm\{load\}\}\}\{\\eta\_\{\\mathrm\{K\}\}\},\\hskip 10\.00002pta\_\{1\}=\-\\frac\{E\}\{\\eta\_\{\\mathrm\{K\}\}\}\.\(25\)For each trial orderα\\alpha, the linear parametersa0a\_\{0\}anda1a\_\{1\}are fitted by least squares, and a bounded one\-dimensional search selectsα\\alpha\. Shape\-preserving interpolation is used only for quadrature; no pointwise fractional derivative is evaluated\. Because the units ofηK\\eta\_\{\\mathrm\{K\}\}depend onα\\alpha, Table[10](https://arxiv.org/html/2608.12879#S5.T10)also reports the dimensionally comparable retardation timetret=\(ηK/E\)1/αt\_\{\\mathrm\{ret\}\}=\(\\eta\_\{\\mathrm\{K\}\}/E\)^\{1/\\alpha\}\. The silt estimate reproduces the reference order and time scale closely\. For clay, the order differs by0\.0980\.098andtrett\_\{\\mathrm\{ret\}\}by about78%78\\%, indicating substantial parameter uncertainty despite a small integral residual\. This section is parameter identification within the prescribed Kelvin support\.

Table 10:Integral\-form identification of the fractional Kelvin model from naturally noisy frozen\-soil creep records\. “Reference” denotes the published constitutive fit\. HereEEis in MPa,ηK\\eta\_\{\\mathrm\{K\}\}has units MPa⋅\\cdottimeα, andtret=\(ηK/E\)1/αt\_\{\\mathrm\{ret\}\}=\(\\eta\_\{\\mathrm\{K\}\}/E\)^\{1/\\alpha\}is in the time unit of the source data; therefore rawηK\\eta\_\{\\mathrm\{K\}\}values at differentα\\alphashould not be compared directly\. The final column is the relative residual of Eq\. \([25](https://arxiv.org/html/2608.12879#S5.E25)\)\.
#### 5\.8Extension to two spatial dimensions

The preceding discovery benchmarks use one spatial coordinate\. To test whether the candidate encoding is tied to that setting, we extend each linear spatial term by a direction labeldj∈\{x,y\}d\_\{j\}\\in\\\{x,y\\\},

Dtα0C​u=∑j=1cξj​𝒳βj\(dj\)​u\.\{\}^\{C\}\_\{0\}D\_\{t\}^\{\\alpha\}u=\\sum\_\{j=1\}^\{c\}\\xi\_\{j\}\\,\\mathcal\{X\}\_\{\\beta\_\{j\}\}^\{\(d\_\{j\}\)\}u\.\(26\)The temporal and two spatial test bases remain separable, and the discrete adjoint is applied along the coordinate named bydjd\_\{j\}\. Mode\-wise tensor contractions avoid assembling a dense spatial Kronecker matrix\. The Pareto search, elbow rule, and exact\-order refit are otherwise unchanged\. The present example restricts the admitted powers top=0p=0and focuses on directional encoding with linear terms\.

We use two doubly periodic benchmarks on\[0,2π\)2\[0,2\\pi\)^\{2\}:

Dt0\.850C​u\\displaystyle\{\}^\{C\}\_\{0\}D\_\{t\}^\{0\.85\}u=0\.30​Dx1\.70​u\+0\.20​Dy1\.40​u,\\displaystyle=0\.30D\_\{x\}^\{1\.70\}u\+0\.20D\_\{y\}^\{1\.40\}u,\(27\)Dt0\.850C​u\\displaystyle\{\}^\{C\}\_\{0\}D\_\{t\}^\{0\.85\}u=−0\.60∂xu\+0\.30Dx1\.70u\+0\.20Dy1\.40u\.\\displaystyle=\-0\.60\\,\\partial\_\{x\}u\+0\.30D\_\{x\}^\{1\.70\}u\+0\.20D\_\{y\}^\{1\.40\}u\.\(28\)Benchmark \([28](https://arxiv.org/html/2608.12879#S5.E28)\) requires the selector to separate two orders in thexxdirection while assigning a third order toyy\. Each Fourier mode is propagated semi\-analytically asu^​\(𝜿,t\)=u^​\(𝜿,0\)​Eα,1​\(λ⁡\(𝜿\)​tα\)\\widehat\{u\}\(\\bm\{\\kappa\},t\)=\\widehat\{u\}\(\\bm\{\\kappa\},0\)E\_\{\\alpha,1\}\(\\lambda\(\\bm\{\\kappa\}\)t^\{\\alpha\}\), independently of the L1 evaluator used by discovery\. On temporal grids withnt=90,179,n\_\{t\}=90,179,and357357, an independent L1 residual is6\.02×10−36\.02\\times 10^\{\-3\},2\.66×10−32\.66\\times 10^\{\-3\}, and1\.18×10−31\.18\\times 10^\{\-3\}, respectively, giving an observed rate1\.181\.18close to the expected2−α=1\.152\-\\alpha=1\.15\.

The reported grid is90×80×8090\\times 80\\times 80\. Applying the paper test\-count rule independently to both spatial axes gives30×40×40=48,00030\\times 40\\times 40=48\{,\}000overlapping weak equations from the same dense data tensor\. Their tensor\-product refinement is computationally tractable because the example reduces repeated objective evaluations to precomputed Gram tables\. We use multiplicative\-uniform noise at0%,1%,5%,10%,0\\%,1\\%,5\\%,10\\%,and20%20\\%,cmax=4c\_\{\\max\}=4,α∈\[0\.55,1\.25\]\\alpha\\in\[0\.55,1\.25\],β∈\[0\.50,2\.50\]\\beta\\in\[0\.50,2\.50\], and the same optimisation and recovery rules as Section[5\.1](https://arxiv.org/html/2608.12879#S5.SS1)\. Boundary trimming is not used\. Table[11](https://arxiv.org/html/2608.12879#S5.T11)reports recovery across all 25 noise–seed combinations\.

Table 11:Two\-dimensional directional discovery over five seeds and five noise levels \(2525runs per benchmark\)\. Support/direction recovery requires the correct support size and direction multiset; complete operator recovery additionally requires all orders to satisfy the protocol of Section[5\.1](https://arxiv.org/html/2608.12879#S5.SS1)\.Both benchmarks retain complete direction and operator recovery throughout the noise sweep\. At20%20\\%noise,\(eα,eβmax,eξmax\)\(e\_\{\\alpha\},e\_\{\\beta\}^\{\\max\},e\_\{\\xi\}^\{\\max\}\)is\(0\.00081±0\.00054,0\.00217±0\.00140,0\.0111±0\.0023\)\(0\.00081\\pm 0\.00054,0\.00217\\pm 0\.00140,0\.0111\\pm 0\.0023\)for Benchmark \([27](https://arxiv.org/html/2608.12879#S5.E27)\) and\(0\.00207±0\.00133,0\.01221±0\.00259,0\.04798±0\.00890\)\(0\.00207\\pm 0\.00133,0\.01221\\pm 0\.00259,0\.04798\\pm 0\.00890\)for Benchmark \([28](https://arxiv.org/html/2608.12879#S5.E28)\)\. A second spatial resolution,90×112×11290\\times 112\\times 112, increases the number of weak rows to94,08094\{,\}080and again gives25/2525/25recoveries for Benchmark \([27](https://arxiv.org/html/2608.12879#S5.E27)\)\. At20%20\\%noise,eξmaxe\_\{\\xi\}^\{\\max\}decreases from0\.01110\.0111to0\.00840\.0084, while the order errors remain of comparable magnitude; we therefore use this result as a grid\-stability check\. Appendix[16](https://arxiv.org/html/2608.12879#S16)reports the resolution and window\-width diagnostics\. The experiment establishes that the encoding and tensor construction extend to coordinate\-dependent orders\.

#### 5\.9Runtime and computational cost

Table[12](https://arxiv.org/html/2608.12879#S5.T12)reports uncached, serial runtimes and search budgets\. The main one\-off cost is precomputing the order\-indexed weak features; each subsequent differential\-evolution evaluation solves a small ridge problem on theKKweak rows\. On the tested benchmarks, both weak and strong frameworks run in a few seconds per seed\. Measurements were obtained in CPU\-only mode with Python 3\.11\.15 on an Apple M4 Pro system with 64 GB unified memory\. They are descriptive implementation\-level timings\. Appendix[17](https://arxiv.org/html/2608.12879#S17)uses the same CPU environment for the neural fractional\-discovery framework, allowing a fair wall\-clock comparison for those particular implementations\.

Table 12:Runtime and search budget on three representative benchmarks \(uncached serial single\-seed runs;KKweak rows; differential\-evolution \(DE\) population multiplier×\\timesgenerations\)\. Times are wall\-clock seconds\.
#### 5\.10Ablation: weak versus strong\-form library, and fixed dictionaries

Table[13](https://arxiv.org/html/2608.12879#S5.T13)separates the three components of Weak\-Pareto: weak measurements, continuous\-order subset search, and exact\-order polishing\. The library\-only contrast is decisive\. With the same continuous\-order selector and no polishing, Strong\-Pareto achieves no complete FADE operator recovery at10%10\\%noise, whereas Weak\-Pareto recovers all five\. With the optimiser held fixed, the0/50/5\-to\-5/55/5change isolates the weak library as the source of the robustness gain\.

A fixed order dictionary is the natural alternative to continuous\-order search, but it faces a resolution–conditioning trade\-off\. A coarse grid introduces order\-discretisation error; a fine grid creates many nearly duplicate columns\. Table[14](https://arxiv.org/html/2608.12879#S5.T14)quantifies this effect for FADE\. AtΔ​β=0\.25\\Delta\\beta=0\.25, the normalised weak dictionary already has condition number2\.6×1092\.6\\times 10^\{9\}; atΔ​β=0\.10\\Delta\\beta=0\.10, its mutual coherence is0\.9870\.987and its condition number reaches8\.4×10158\.4\\times 10^\{15\}\. Such a large condition number indicates severe numerical ill\-conditioning, or near rank deficiency: small perturbations in the data or arithmetic can cause large changes in fitted coefficients\. It is therefore a numerical indicator of regression instability\. Weak\-Pareto avoids this dense design by evaluating at mostc≤cmaxc\\leq c\_\{\\max\}proposed columns at a time\. The final direct evaluation removes interpolation error conditional on the selected model, while support selection can still depend on the order grid and search trajectory\.

The empirical ablation confirms the theoretical distinction\. Weak Grid\-STRidge, which applies sequential threshold ridge regression \(STRidge\) on a fixed weak dictionary, does not recover the correct FADE support in any of the five runs even though its fixed grid contains nodes within0\.010\.01of both true orders\. Its small residual therefore reflects a coherent overcomplete dictionary that fits the weak equations without identifying the correct structure\. Continuous\-order Weak\-Pareto recovers all five supports and operators\. Disabling exact\-order polishing leaves recovery unchanged and changeseβmaxe\_\{\\beta\}^\{\\max\}only from0\.080\.08to0\.070\.07; polishing is therefore a small refinement on this benchmark\. Together, the ablations show that weak measurements provide noise robustness and continuous\-order best\-subset search provides reliable support selection\.

Table 13:Component ablation on FADE at10%10\\%noise, toggling the weak library, continuous\-order search, and exact\-order polishing\.Table 14:Conditioning of a fixed dense fractional dictionary as the order spacingΔ​β\\Delta\\betashrinks: number of columns, mutual coherence, and condition number of the column\-normalised weak library \(noiseless FADE, orders over\[0\.1,3\]\[0\.1,3\]\)\.
#### 5\.11Limitations and future work

Two empirical limits are prominent\. In the Riesz reaction–diffusion experiments, support and powers are often retained while nearby spatial orders and small reaction coefficients remain difficult to distinguish\. Increasing the optimisation budget does not resolve this ambiguity, and Section[3\.5](https://arxiv.org/html/2608.12879#S3.SS5)shows that these spatial orders have substantially flatter coefficient\-profiled directions than the directional FADE and Burgers terms\. Their accurate clean profiles indicate sensitivity to perturbations once noise is present\. In the fixed\-support superunit diagnostic, the correct branch is retained at1%1\\%noise, but four of five temporal\-order estimates reach the upper search bound\. Remark[2](https://arxiv.org/html/2608.12879#Thmremark2)links this behaviour to the endpoint\-concentrated composed L1 target; alternative treatments of the initial rate are left for future work\.

The theoretical guarantees are strongest for linear right\-hand\-side features\. Nonlinear weak features average the corresponding strong features but reuse the noisy field and can therefore be biased\. The support\-conditioned sensitivity result is local; global identifiability of the joint discrete–continuous model class remains open\. The differential\-evolution search and local polishing are heuristic optimisation procedures, and rigorous global\-convergence analysis is left for future work\.

Several empirical directions remain open\. Test\-window scale is problem dependent: in the one\-dimensionalKK\-sweep, increasing the number of rows narrows the localised Gaussian test\-function windows and eventually degrades operator recovery, while substantially broader Gaussian windows cause support under\-selection in the two\-dimensional diagnostic\. The challenging Riesz cases also depend on the admitted powers and selection rule\. The two\-dimensional study establishes directional extensibility for dense, periodic, uniformly sampled data with48,00048\{,\}000overlapping weak rows\. Sparse observations, nonlinear multidimensional discovery, and matched two\-dimensional weak–strong comparisons remain future directions\. Observation horizon and spectral content were not varied systematically, and overlapping test functions leave the disjoint validation rows statistically correlated\.

Accordingly, we report support recovery, operator recovery, order error, coefficient error, and fit residual separately\. Operator\-specific care also remains essential: finite\-domain and periodic fractional derivatives are different models, and each weak feature must use the adjoint of the operator it represents\.

### 6Conclusion

Weak\-Pareto combines two methodological advances for fractional equation discovery: an adjoint\-consistent weak library and a continuous\-order Pareto search\. For linear right\-hand\-side terms, the weak formulation removes pointwise fractional differentiation from the measured field, and the variance analysis explains why this improves robustness as the grid is refined\. For nonlinear terms, it provides integrated projections that reduce noise through averaging but may remain biased\. The continuous\-order encoding avoids the discretisation error and severe collinearity associated with fixed dictionaries, while the elbow search makes the trade\-off between fit and complexity explicit\. A support\-conditioned local sensitivity analysis additionally quantifies when spatial\-order changes can be absorbed by coefficient refitting, linking weak regression geometry to the observed hierarchy of noisy spatial\-order errors\.

The experiments support these advantages\. On FADE and fractional Burgers, Weak\-Pareto recovers the correct support in every seed at all tested multiplicative\-noise levels up to20%20\\%\. The matched strong\-form framework largely fails once noise is introduced, and the conclusion is unchanged under additive Gaussian noise\. Component ablations show that the weak library drives noise robustness and that continuous\-order search resolves the support\-selection failure of the fixed\-grid baseline on FADE\. The superunit diagnostic gives 5/5 operator recovery through0\.5%0\.5\\%noise but 1/5 at1%1\\%, marking an endpoint\-sensitive limit\. On the advection–diffusion benchmark, Weak\-Pareto also yields more consistent operator recovery and lower wall\-clock runtime than the adapted neural fractional\-discovery framework under the same CPU environment, although the operator realisations are not identical\. Together, the Riesz and superunit results show that support or branch selection can remain stable while fractional orders become weakly identifiable\. The two\-dimensional example further shows that the candidate tuple can absorb coordinate direction as an additional discrete label: both anisotropic benchmarks are recovered in every run through20%20\\%noise\. Finally, the frozen\-soil example shows that the weak\-form representation can fit a fractional Kelvin model to irregular experimental data, with close agreement for silt and greater uncertainty for clay\. Future work should address uncertainty\-aware selection near identifiability limits, sparse and nonlinear multidimensional discovery, and application to experimental anomalous\-transport systems\.

### 7Adjoint identities

##### Riemann–Liouville integration by parts\.

For0<γ<10<\\gamma<1and sufficiently regularf,ϕf,\\phi, the left RL derivative \([1](https://arxiv.org/html/2608.12879#S2.E1)\) satisfies

∫ab\(Dγza​f\)​\(z\)​ϕ​\(z\)​𝑑z=∫abf⁡\(z\)​\(Dγbz​ϕ\)​\(z\)​𝑑z,\\int\_\{a\}^\{b\}\(\{\}\_\{a\}D\_\{z\}^\{\\gamma\}f\)\(z\)\\,\\phi\(z\)\\,\\mathrm\{d\}z=\\int\_\{a\}^\{b\}f\(z\)\\,\(\{\}\_\{z\}D\_\{b\}^\{\\gamma\}\\phi\)\(z\)\\,\\mathrm\{d\}z,\(29\)i\.e\. a left derivative onffbecomes a right derivative onϕ\\phi\. For0<γ<10<\\gamma<1it suffices in the present application thatffbe bounded andϕ\\phivanish at the right endpointbb\. More generally, iff∈Lq​\(a,b\)f\\in L^\{q\}\(a,b\)withq\>1/\(1−γ\)q\>1/\(1\-\\gamma\), then\(I1−γza​f\)​\(z\)=O⁡\(\(z−a\)1−γ−1/q\)→0\(\{\}\_\{a\}I\_\{z\}^\{1\-\\gamma\}f\)\(z\)=O\(\(z\-a\)^\{1\-\\gamma\-1/q\}\)\\to 0asz→a\+z\\to a^\{\+\}; boundedffis a sufficient special case and applies here tof=u−u⁡\(0,⋅\)f=u\-u\(0,\\cdot\)\. Withϕ⁡\(b\)=0\\phi\(b\)=0, one hasdd​z​\(I1−γbz​ϕ\)=I1−γbz​\(ϕ′\)\\frac\{\\,\\mathrm\{d\}\}\{\\,\\mathrm\{d\}z\}\\bigl\(\{\}\_\{z\}I\_\{b\}^\{1\-\\gamma\}\\phi\\bigr\)=\{\}\_\{z\}I\_\{b\}^\{1\-\\gamma\}\(\\phi^\{\\prime\}\), from which \([29](https://arxiv.org/html/2608.12879#S7.E29)\) follows by ordinary integration by parts and Fubini\.

For higher ordersn−1<γ<nn\-1<\\gamma<n, the analogous identity also requires the relevant left\-endpoint traces ofIn−γza​f\{\}\_\{a\}I\_\{z\}^\{n\-\\gamma\}fto vanish andϕ\(m\)​\(b\)=0\\phi^\{\(m\)\}\(b\)=0form=0,…,n−1m=0,\\dots,n\-1\. In the present superunit Caputo application,f=u−P1,a​uf=u\-P\_\{1,a\}usatisfiesf⁡\(a\)=f′​\(a\)=0f\(a\)=f^\{\\prime\}\(a\)=0; for example,u⁡\(⋅,x\)∈C1​\[a,b\]u\(\\cdot,x\)\\in C^\{1\}\[a,b\]with locally Hölder\-continuousutu\_\{t\}is a sufficient regularity condition for these left traces to vanish\.

If an endpoint condition fails, its continuum boundary term must be retained\. The reported Gaussian tests instead use the exact discrete\-adjoint construction of Section[3\.3](https://arxiv.org/html/2608.12879#S3.SS3), which preserves the implemented discrete inner\-product identity without assuming vanishing Gaussian traces\.

##### Caputo target\.

WritingDtα0C​u=Dαt0​\[u−Pn−1,0​u\]\{\}\_\{0\}^\{C\}\\\!D\_\{t\}^\{\\alpha\}u=\{\}\_\{0\}D\_\{t\}^\{\\alpha\}\[u\-P\_\{n\-1,0\}u\]and applying the higher\-order counterpart of Eq\. \([29](https://arxiv.org/html/2608.12879#S7.E29)\) gives Eq\. \([9](https://arxiv.org/html/2608.12879#S3.E9)\) for0<α<20<\\alpha<2,α≠1\\alpha\\neq 1, provided the test function satisfies the corresponding terminal conditions\. For0<α<10<\\alpha<1, the subtracted polynomial isu⁡\(0,⋅\)u\(0,\\cdot\); for1<α<21<\\alpha<2, it isu⁡\(0,⋅\)\+t​∂tu⁡\(0,⋅\)u\(0,\\cdot\)\+t\\,\\partial\_\{t\}u\(0,\\cdot\)\. The discrete target implements the same branch\-specific correction through the transpose of the corresponding Caputo matrix\.

##### Spectral operators\.

For periodic fields, Parseval’s identity and the \(conjugate\) symmetry of the multipliers in \([2](https://arxiv.org/html/2608.12879#S2.E2)\) give self\-adjointness of the Riesz operator and the conjugate\-multiplier adjoint of the directional operator stated in Section[3\.2](https://arxiv.org/html/2608.12879#S3.SS2)\.

### 8Discrete operators and numerical verification

The discrete adjoint identity \([12](https://arxiv.org/html/2608.12879#S3.E12)\) is verified for every operator family: matrix transposes for Grünwald–Letnikov and one\-sided finite\-domain stencils, the transposed L1 matrix for Caputo time derivatives, and conjugate Fourier multipliers for periodic Riesz and directional operators\. Random smooth\-field tests satisfy⟨A​f,ϕ⟩h=⟨f,A∗,h​ϕ⟩h\\langle Af,\\,\\phi\\rangle\_\{h\}=\\langle f,\\,A^\{\\ast,h\}\\phi\\rangle\_\{h\}to machine precision\. The optional fractional\-integral adjoint is verified in the same way\.

We also checked the grid\-refinement assumptions of Proposition[1](https://arxiv.org/html/2608.12879#Thmtheorem1)empirically\. Fornt=nx=n∈\{24,32,48,64,96\}n\_\{t\}=n\_\{x\}=n\\in\\\{24,32,48,64,96\\\}, a fixed smooth periodic separable test function, a spectral orderβ=1\.5\\beta=1\.5, and20002000independent standard\-Gaussian noise fields per grid, least\-squares fits of log variance againstlog⁡n\\log ngave slopes−2\.01\-2\.01for the weak feature and2\.972\.97for the pointwise feature, close to the predicted−2\-2and2​β=32\\beta=3\. An independent deterministic calculation reproduces this check\. This verification holds the test function onQQfixed, exactly as assumed in the proposition; the separateKK\-sensitivity study of Appendix[12](https://arxiv.org/html/2608.12879#S12)examines what happens when the number and width of the localised Gaussian test functions are changed\.

The exact adjoint check verifies the discrete transpose identities\. Continuum consistency is assessed independently using analytic operator–function pairs\. Forℛβ​sin⁡\(m​x\)=−\|m\|β​sin⁡\(m​x\)\\mathcal\{R\}\_\{\\beta\}\\sin\(mx\)=\-\|m\|^\{\\beta\}\\sin\(mx\)withβ=1\.7\\beta=1\.7andm=3m=3, the periodic FFT implementation has relative errors between8×10−158\\times 10^\{\-15\}and2×10−132\\times 10^\{\-13\}on grids of3232–256256points\. ForDt0\.70C​t3=Γ⁡\(4\)​t2\.3/Γ⁡\(3\.3\)\{\}\_\{0\}^\{C\}D\_\{t\}^\{0\.7\}t^\{3\}=\\Gamma\(4\)t^\{2\.3\}/\\Gamma\(3\.3\), the Caputo–L1 relative error decreases from5\.2×10−35\.2\\times 10^\{\-3\}to3\.6×10−43\.6\\times 10^\{\-4\}over6565–513513points, with observed rates1\.281\.28–1\.291\.29, close to the expected2−α=1\.32\-\\alpha=1\.3\. For the separately implemented superunit composition,Dt1\.30C​t3=Γ⁡\(4\)​t1\.7/Γ⁡\(2\.7\)\{\}\_\{0\}^\{C\}D\_\{t\}^\{1\.3\}t^\{3\}=\\Gamma\(4\)t^\{1\.7\}/\\Gamma\(2\.7\), the relative error over the non\-initial rows decreases from1\.20×10−21\.20\\times 10^\{\-2\}to9\.99×10−49\.99\\times 10^\{\-4\}over the same grids, with observed rates1\.191\.19–1\.201\.20\. We claim empirical convergence for this composed discretisation; a formal convergence order for the complete composition is not derived here\. This verifies convergence of the active superunit code path independently of model selection\. The semi\-analytic fixed\-support recovery experiment in Section[5\.3](https://arxiv.org/html/2608.12879#S5.SS3)provides the complementary branch/order test on clean and noisy fields\.

The nonlinear\-bias formulas of Remark[1](https://arxiv.org/html/2608.12879#Thmremark1)were also checked numerically\. For the periodic first derivative, the predicted leading bias is zero\. The observed Monte Carlo means remain close to zero and, unlike the order\-1\.71\.7results, show no systematic growth with grid refinement\. For a directional derivative of order1\.71\.7, the additive\-Gaussian predicted biases atnx=64,128,256n\_\{x\}=64,128,256are−8\.39\-8\.39,−27\.25\-27\.25, and−88\.52\-88\.52, while the observed means are−8\.36\-8\.36,−27\.23\-27\.23, and−88\.60\-88\.60\. Under multiplicative\-uniform noise the corresponding predictions are−1\.85\-1\.85,−6\.02\-6\.02, and−19\.55\-19\.55, versus observed means−1\.85\-1\.85,−6\.03\-6\.03, and−19\.58\-19\.58\. The successive additive\-bias ratios are about3\.253\.25, matching the predicted grid\-doubling factor21\.72^\{1\.7\}\. The observed growth agrees with the predicted nonlinear\-feature bias, confirming an intrinsic limitation of these data\-weighted nonlinear features\.

Finally, an oracle diagnostic profiles the held\-out objective on the correct two\-term support and operator modes while minimising over the positive Riesz order\. Its five arms use clean data, fully noisy data, the noisy field with only the clean initial slice restored, a noisy target with a clean library, and a clean target with a noisy library\. For the space\-fractional benchmark the branch\-aware clean profile selects the exact integer modeα=1\\alpha=1; at2%2\\%noise the fully noisy and target\-only arms minimise at0\.8200\.820and0\.8170\.817, whereas the library\-only arm remains at the exact integer mode\. Restoring only the initial slice gives0\.8000\.800, the declared lower bound; the constrained profile therefore remains shifted and does not provide an interior minimum\. For the time–space benchmark the clean and library\-only minima are0\.8200\.820, while the fully noisy and target\-only minima are both0\.7720\.772; restoring the initial slice gives0\.7630\.763\. This fixed\-support diagnostic identifies the noisy temporal target as the dominant source of the observed shift in these examples\. The direction of that shift is profile\-specific and need not persist for other problems or scoring conventions\.

### 9Experimental settings

The one\-dimensional benchmark grids are listed in Table[2](https://arxiv.org/html/2608.12879#S5.T2); all spatial grids are periodic, and their sizes range from80×8080\\times 80to150×120150\\times 120in time–space samples\. The two\-dimensional example uses a90×80×8090\\times 80\\times 80time–space–space grid\. The reported weak library uses tensor\-product Gaussian test functions for most searches and Fourier spatial modes for the periodic high\-order Riesz cases\. The main configuration uses powersp∈\{0,1,2\}p\\in\\\{0,1,2\\\},cmax=4c\_\{\\max\}=4, automatic plateau and selection\-stability stopping, a relative\-improvement toleranceδ=0\.03\\delta=0\.03, and the signed elbow rule with two\-point margin0\.150\.15\. Appendix[12](https://arxiv.org/html/2608.12879#S12)examines sensitivity to the number of weak rows\.

Fractional orders are represented on uniform, branch\-confined grids: 47 nominal temporal nodes and 59 spatial nodes over the declared ranges in Table[15](https://arxiv.org/html/2608.12879#S9.T15)\. The separatorϵα=10−3\\epsilon\_\{\\alpha\}=10^\{\-3\}prevents interpolation acrossα=1\\alpha=1, which is represented by a distinct exact\-integer candidate\. No noninteger true order is inserted into the grids or optimiser initialisation\. The identity is also a distinct spatial candidate atβ=0\\beta=0\. After selection, direct operator evaluation at the polished orders removes interpolation error from the conditional refit, but the selected model can still depend on grid resolution and the search trajectory\.

Differential evolution uses a SciPy population multiplier of77and 24 generations\. The ridge parameter isλ=10−3\\lambda=10^\{\-3\}on column\-normalised designs, the validation fraction is0\.250\.25, and pruning usesτcontrib=10−4\\tau\_\{\\mathrm\{contrib\}\}=10^\{\-4\}andτabs=10−10\\tau\_\{\\mathrm\{abs\}\}=10^\{\-10\}, with a relative\-coefficient fallback of10−310^\{\-3\}when library columns are unavailable\. The residual guards are10−1410^\{\-14\}for validation scoring and10−1210^\{\-12\}for coefficient errors\. Temporal boundary trimming is a control for pointwise feature construction only; it is not applied to the weak framework, whose discrete adjoint target retains the endpoint structure of the implemented operator\. For the nonnegative FADE and integer advection–diffusion \(ADE\) fields, incorrect nonlinear candidate terms use\(u\+\)p\(u\_\{\+\}\)^\{p\}to avoid amplifying noise\-induced sign changes; the sign\-changing Burgers field usesupu^\{p\}\. The challenging cases in Appendix[11](https://arxiv.org/html/2608.12879#S11)use the stated case\-specific power set and, for the two\-term Riesz case, the heuristic Akaike information criterion \(AIC\)\-type selector\.

Table 15:Declared search domains per benchmark: the temporal\-order interval\[αmin,αmax\]\[\\alpha\_\{\\min\},\\alpha\_\{\\max\}\], the spatial\-order interval\[βmin,βmax\]\[\\beta\_\{\\min\},\\beta\_\{\\max\}\], the support\-size capcmaxc\_\{\\max\}, and the admitted power set𝒜\\mathcal\{A\}\. These are the differential\-evolution bounds \(Section[4](https://arxiv.org/html/2608.12879#S4)\); they encode coarse prior knowledge of the operator regime\. No noninteger benchmark\-true fractional order is inserted, whereas exact integer modes are included independently of the benchmark truth\. If a temporal interval crosses one, its fractional portions and the exact integer candidate are searched separately\. The identity operator is a discrete candidate atβ=0\\beta=0; derivative orders are searched on the strictly positive part of the interval\.For ADE, the classical second derivativeβ=2\\beta=2is simply the upper integer endpoint of the search range\. A supplementary reduced\-sampling diagnostic retains every second temporal snapshot without interpolation or imputation and is included only as an implementation check\.

##### Local spatial\-order diagnostic\.

Eq\. \([16](https://arxiv.org/html/2608.12879#S3.E16)\) is evaluated with a centred finite difference of step10−410^\{\-4\}on exact weak features\. The fixed\-support profiles use 61 equally spaced trial orders over the benchmark’s positive spatial\-order search interval, augmented by the true order, followed by bounded one\-dimensional refinement around the best grid point with tolerance10−510^\{\-5\}\. The true temporal branch/order, support, and all other spatial orders are fixed; coefficients, training/validation rows, ridge parameter, and validation score follow the main protocol\. The reproduction script and all0%0\\%,2%2\\%,5%5\\%, and10%10\\%per\-seed profiles are archived in Online Resource 1\.

#### 9\.1Implementation of the best\-subset selection

Eqs\. \([19](https://arxiv.org/html/2608.12879#S4.E19)\)–\([20](https://arxiv.org/html/2608.12879#S4.E20)\) define the outer objective and inner coefficient fit\. Weak rows are split deterministically into training and validation subsets, with validation fraction0\.250\.25\. Training columns are normalised before ridge regression and the fitted coefficients are mapped back to the original scale\. The principal criterion for differential evolution, Pareto dominance, elbow selection, and stopping is the variance\-normalised validation score in Eq\. \([21](https://arxiv.org/html/2608.12879#S4.E21)\)\. Training error, validation error, the selection objective, and the post\-refit full\-data residual remain distinct fields: the objective is normallylog10\\log\_\{10\}of the normalised validation mean\-squared error \(MSE\) plus any declared search penalty, whereasℰfit\\mathcal\{E\}\_\{\\mathrm\{fit\}\}is computed only after exact\-order refitting\. The duplicate\-order penalty is

Πdup=λdup​∑i<jpi=pj\(1−\|βi−βj\|δβ\)\+,λdup=0\.02,δβ=0\.04\.\\Pi\_\{\\mathrm\{dup\}\}=\\lambda\_\{\\mathrm\{dup\}\}\\\!\\sum\_\{\\begin\{subarray\}\{c\}i<j\\\\ p\_\{i\}=p\_\{j\}\\end\{subarray\}\}\\\!\\left\(1\-\\frac\{\|\\beta\_\{i\}\-\\beta\_\{j\}\|\}\{\\delta\_\{\\beta\}\}\\right\)\_\{\+\},\\hskip 20\.00003pt\\lambda\_\{\\mathrm\{dup\}\}=0\.02,\\hskip 10\.00002pt\\delta\_\{\\beta\}=0\.04\.\(30\)The branch separator isϵα=10−3\\epsilon\_\{\\alpha\}=10^\{\-3\}\. Differential evolution is run independently on each nonempty temporal mode, using the branch\-specific dimensions and bounds stated in Section[4\.2](https://arxiv.org/html/2608.12879#S4.SS2); its best candidates are then compared by the same objective\. The exact integer mode evaluates∂t\\partial\_\{t\}directly and never interpolates neighbouring fractional features\. The selected mode is retained during inactive\-term pruning, exact\-order polishing, and the final full\-row coefficient refit\.

### 10Nonlinear fractional Burgers solver

The nonlinear benchmark∂tu=−u∂xu\+νDxβu\\partial\_\{t\}u=\-u\\,\\partial\_\{x\}u\+\\nu D\_\{x\}^\{\\beta\}uis integrated pseudospectrally on480480periodic spatial points over\[0,30\)\[0,30\), using fourth\-order Runge–Kutta with internal stepΔ​tfine=0\.004\\Delta t\_\{\\mathrm\{fine\}\}=0\.004up toT=12T=12\. The quadratic flux−12∂x\(u2\)\-\\tfrac\{1\}\{2\}\\partial\_\{x\}\(u^\{2\}\)is dealiased by the2/32/3rule\. Retaining every fourth spatial point and150150endpoint\-excluded temporal snapshots gives the reported150×120150\\times 120grid\. The conjugate\-symmetric multiplier keeps the field real, whileRe⁡\(i​κ\)β<0\\operatorname\{Re\}\(\\mathrm\{i\}\\kappa\)^\{\\beta\}<0for1<β<21<\\beta<2supplies dissipation\. The parametersν=0\.25\\nu=0\.25andβ=1\.7\\beta=1\.7keep the solution smooth while maintaining comparable nonlinear and diffusive contributions\.

### 11Additional experiments: challenging cases and hyperparameter sensitivity

We report two cases that probe the limits of the method under the non\-restrictive main\-text settings \(plateau stopping on,cmax=4c\_\{\\max\}=4, powersp∈\{0,1,2\}p\\in\\\{0,1,2\\\}, elbow selection\)\. They are at opposite ends of the model class:

- •No fractional derivative \(ADE\)\.The integer\-order advection–diffusion equation∂tu=−∂xu\+0\.25∂x2u\\partial\_\{t\}u=\-\\partial\_\{x\}u\+0\.25\\,\\partial\_\{x\}^\{2\}u, with true terms\(0,1,−1\)\(0,1,\-1\)and\(0,2,0\.25\)\(0,2,0\.25\)\(directional operator\)\.
- •More than one fractional spatial derivative \(two\-term Riesz\)\.The equation∂tu=0\.05​ℛ0\.55​u\+0\.005​ℛ2\.8​u\\partial\_\{t\}u=0\.05\\,\\mathcal\{R\}\_\{0\.55\}u\+0\.005\\,\\mathcal\{R\}\_\{2\.8\}u, with true terms\(0,0\.55,0\.05\)\(0,0\.55,0\.05\)and\(0,2\.8,0\.005\)\(0,2\.8,0\.005\)\.

Under the common main\-text settings, both challenging cases fail structurally \(Table[16](https://arxiv.org/html/2608.12879#S11.T16)\)\. For ADE, the search selects the correct support size but replaces∂x2u\\partial\_\{x\}^\{2\}uwith a nonlinear candidate that is absent from the true equation\. For the two\-term Riesz equation, the dominant high\-order term is identified, but the low\-order term contributes less than the noise floor at10%10\\%; consequently, the elbow selects only one term\. Because the support is then incomplete, order and coefficient errors are not reported for the default rows\.

Case\-specific restrictions improve support recovery but not complete operator identification\. Limiting ADE top=0p=0recovers the support and powers in4/54/5seeds\. For the two\-term Riesz case, combiningp=0p=0with a heuristic AIC\-type selector recovers the two\-term support in all five seeds\. This selector usesnval​log⁡MSEval\+2​kn\_\{\\mathrm\{val\}\}\\log\\mathrm\{MSE\}\_\{\\mathrm\{val\}\}\+2kon held\-out rows\. Because changingα\\alphachanges the response, this score is used only as a case\-specific heuristic; the variance\-normalised validation score remains the principal selector\. Neither adjustment achieves complete operator recovery at10%10\\%noise, confirming sensitivity to the admitted powers, selection rule, and order identifiability\.

Table 16:Challenging cases under non\-restrictive defaults versus case\-specific adjustments at10%10\\%noise\. Errors are conditioned on support/power recovery;ℰfit\\mathcal\{E\}\_\{\\mathrm\{fit\}\}is over all seeds\.#### 11\.1Two\-point elbow\-margin sensitivity

The special two\-point rule in Section[4\.3](https://arxiv.org/html/2608.12879#S4.SS3)uses the marginm2=0\.15m\_\{2\}=0\.15only when the available front containsc=1c=1andc=2c=2\. Define

Δ12=log10⁡ℰ1⋆−log10⁡ℰ2⋆=log10⁡\(ℰ1⋆/ℰ2⋆\),\\Delta\_\{12\}=\\log\_\{10\}\\mathcal\{E\}\_\{1\}^\{\\star\}\-\\log\_\{10\}\\mathcal\{E\}\_\{2\}^\{\\star\}=\\log\_\{10\}\(\\mathcal\{E\}\_\{1\}^\{\\star\}/\\mathcal\{E\}\_\{2\}^\{\\star\}\),so thatc=2c=2clears the two\-point rule whenΔ12\>m2\\Delta\_\{12\}\>m\_\{2\}\. Table[17](https://arxiv.org/html/2608.12879#S11.T17)appliesm2∈\{0\.10,0\.15,0\.20\}m\_\{2\}\\in\\\{0\.10,0\.15,0\.20\\\}to the same paper\-budgetc=1,2c=1,2fronts at10%10\\%noise, so the comparison isolates the margin itself\. Thec=2c=2candidates are selection\-stage fits from these restricted two\-point fronts\. Section[5\.2](https://arxiv.org/html/2608.12879#S5.SS2)reports the fully searched and refined models, so the corresponding parameter estimates can differ\.

Table 17:Sensitivity of the two\-point elbow decision at10%10\\%multiplicative noise \(five seeds\)\. The interval gives the observed range ofΔ12\\Delta\_\{12\}across seeds; the final three columns count seeds in whichc=2c=2clears the stated margin\. The reported default ism2=0\.15m\_\{2\}=0\.15\.The FADE, fractional Burgers, and time–space reaction–diffusion two\-point decisions are unchanged throughout this neighbourhood of the default margin\. The space\-fractional Riesz case is selector\-sensitive at10%10\\%noise: full\-selector reruns give support/power recovery of5/55/5,3/53/5, and1/51/5at margins0\.100\.10,0\.150\.15, and0\.200\.20, respectively, while complete operator recovery remains0/50/5throughout\. The margin therefore affects support retention in this severe case, while the spatial\-order difficulty persists\. Online Resource 1 archives the completec=1,2c=1,2audit and the full\-selector verification\.

### 12Sensitivity to the number of weak rows

The number of weak rowsK=nttest​nxtestK=n\_\{t\}^\{\\mathrm\{test\}\}n\_\{x\}^\{\\mathrm\{test\}\}controls measurement resolution, not the candidate class\. Table[18](https://arxiv.org/html/2608.12879#S12.T18)variesKKseventeen\-fold on FADE at10%10\\%noise\. Support and powers are recovered in all five seeds throughout, and complete operator recovery remains5/55/5up to the main settingK=2640K=2640\. AtK=5270K=5270, support recovery remains perfect but operator recovery falls to1/51/5, with larger spatial\-order and coefficient errors\.

At largeKK, the localised Gaussian test\-function windows generated by the paper’s count\-to\-width rule become narrower in time and space; narrower localisation broadens their spectra, allowing more high\-wavenumber noise to enter each weak row\. ModerateKKtherefore improves coefficient precision and reduces construction cost\. The few\-column continuous\-order design remains well conditioned, with condition numberO⁡\(1\)O\(1\), unlike the dense fixed dictionary in Table[14](https://arxiv.org/html/2608.12879#S5.T14)\. The optimal test\-function family and bandwidth remain problem dependent; a matched comparison of Gaussian, compact\-bump, and Fourier families is left to future work\.

Table 18:Sensitivity to the number of weak rowsK=nttest×nxtestK=n\_\{t\}^\{\\mathrm\{test\}\}\\times n\_\{x\}^\{\\mathrm\{test\}\}on FADE at10%10\\%noise \(five seeds\)\. Reported dispersions are sample standard deviations\. Under the Gaussian count\-to\-width rule, the number and width of localised test functions are coupled; largerKKtherefore means narrower Gaussian test\-function windows\. TheK=2640K=2640row is the main setting\.
### 13Alternative noise law

The main experiments use multiplicative uniform perturbations because they preserve the local signal scale\. To test whether the weak\-versus\-strong gap depends on that particular law, Table[19](https://arxiv.org/html/2608.12879#S13.T19)repeats the10%10\\%FADE and fractional Burgers comparisons with independent additive Gaussian noise,

u~=u\+0\.10​std⁡\(u\)​Z,Zi​j∼𝒩⁡\(0,1\),\\widetilde\{u\}=u\+0\.10\\,\\operatorname\{std\}\(u\)Z,\\hskip 20\.00003ptZ\_\{ij\}\\sim\\mathcal\{N\}\(0,1\),using the same five seeds, search ranges, and optimisation budgets as the main experiments\. The same noisy field is supplied to both methods for every seed\. Weak\-Pareto recovers the correct support in all five runs for both equations and the complete operator in4/54/5seeds for each; the strong\-form framework does not recover the correct support in any run for either equation\. The result is consistent with Corollary[2](https://arxiv.org/html/2608.12879#Thmtheorem2): the averaging advantage of the weak library persists under additive Gaussian noise\.

Table 19:Additive\-Gaussian robustness with noise standard deviation equal to10%10\\%of the clean\-field standard deviation \(five seeds\)\. Parameter errors are conditioned on support/power recovery;ℰfit\\mathcal\{E\}\_\{\\mathrm\{fit\}\}is over all seeds\.
### 14Forward\-model validation

Weak residual and trajectory reproduction measure different properties\. We therefore integrate each representative discovered FPDE from the benchmark initial condition and report

efield=∥udisc−u∥2∥u∥2\+ε\.e\_\{\\mathrm\{field\}\}=\\frac\{\\lVert u\_\{\\mathrm\{disc\}\}\-u\\rVert\_\{2\}\}\{\\lVert u\\rVert\_\{2\}\+\\varepsilon\}\.\(31\)The spatial integrator uses the benchmark’s declared periodic Riesz or directional multiplier\. Time integration uses the exact exponential propagator for linear integer\-time equations, an adaptively substepped and dealiased Runge–Kutta scheme for Burgers, and the Caputo L1 scheme for fractional time\. The same solver is also run with the true parameters to quantify numerical discrepancy \(Table[20](https://arxiv.org/html/2608.12879#S14.T20)\)\. For time–space reaction–diffusion, the generator and evaluator share the L1 discretisation, so this row is a self\-consistency check\. FADE is generated semi\-analytically; both discovered\-model errors are no larger than the true\-parameter solver discrepancy, making those rows inconclusive\. All simulations start from the discovery initial condition and cover the same horizon, so the test measures trajectory reproduction\. The ADE control is omitted because the generic periodic evaluator uses a different operator convention from that finite\-domain dataset; even the true\-parameter simulation exceeds the self\-consistency threshold\.

Table 20:Forward\-model validation: normalised field errorefielde\_\{\\mathrm\{field\}\}between the discovered\-model simulation and the clean reference, with the true\-parameter solver discrepancy for context\. This is a same\-initial\-condition trajectory\-reproduction test; held\-out prediction is outside its scope\.
### 15Closed\-form discovered equations

Table[21](https://arxiv.org/html/2608.12879#S15.T21)lists representative clean and noisy discoveries, including the selected powers, orders, and post\-pruning coefficients\. Each representative seed has the median worst spatial\-order error, avoiding a best\-case presentation; a dagger marks failure of the operator\-structure criterion\. The noisy reaction–diffusion representatives attain their lower temporal bounds \(α^=0\.8000\\widehat\{\\alpha\}=0\.8000and0\.65000\.6500\); these boundary\-attaining estimates are consistent with Appendix[8](https://arxiv.org/html/2608.12879#S8)and are not interior optima\.

Table 21:Representative discovered equations versus ground truth\. The representative seed has the medianeβmaxe\_\{\\beta\}^\{\\max\}; a dagger marks an unrecovered operator structure\. Exact integer temporal modes are printed as∂t\\partial\_\{t\}\. Exact integer orders in the ground\-truth equations use derivative shorthand, whereas continuously estimated spatial orders are printed numerically, including estimates that round to1\.001\.00\.BenchmarkNoiseDiscovered equation \(representative seed\)*FADE*— true:Dt0\.80​u=−1\.00​ux\+0\.50​Dx1\.70​uD\_\{t\}^\{0\.80\}u=\-1\.00\\,u\_\{x\}\+0\.50\\,D\_\{x\}^\{1\.70\}uFADE0%Dt0\.7990​u=−0\.99​Dx1\.00​u\+0\.49​Dx1\.73​uD\_\{t\}^\{0\.7990\}u=\-0\.99\\,D\_\{x\}^\{1\.00\}u\+0\.49\\,D\_\{x\}^\{1\.73\}u10%Dt0\.8038​u=−1\.11​Dx1\.03​u\+0\.63​Dx1\.61​uD\_\{t\}^\{0\.8038\}u=\-1\.11\\,D\_\{x\}^\{1\.03\}u\+0\.63\\,D\_\{x\}^\{1\.61\}u*Frac\. RD \(space\)*— true:∂tu=0\.04​u\+0\.18​ℛ1\.65​u\\partial\_\{t\}u=0\.04\\,u\+0\.18\\,\\mathcal\{R\}\_\{1\.65\}uFrac\. RD\(space\)0%∂tu=0\.04​u\+0\.18​ℛ1\.65​u\\partial\_\{t\}u=0\.04\\,u\+0\.18\\,\\mathcal\{R\}\_\{1\.65\}u5%Dt0\.8000​u=−0\.06​ℛ0\.27​u\+0\.20​ℛ1\.53​uD\_\{t\}^\{0\.8000\}u=\-0\.06\\,\\mathcal\{R\}\_\{0\.27\}u\+0\.20\\,\\mathcal\{R\}\_\{1\.53\}u†\\dagger*Frac\. RD \(time–space\)*— true:Dt0\.82​u=0\.03​u\+0\.12​ℛ1\.55​uD\_\{t\}^\{0\.82\}u=0\.03\\,u\+0\.12\\,\\mathcal\{R\}\_\{1\.55\}uFrac\. RD\(time–space\)0%Dt0\.8208​u=0\.03​u\+0\.12​ℛ1\.55​uD\_\{t\}^\{0\.8208\}u=0\.03\\,u\+0\.12\\,\\mathcal\{R\}\_\{1\.55\}u5%Dt0\.6500​u=−0\.06​ℛ0\.29​u\+0\.16​ℛ1\.44​uD\_\{t\}^\{0\.6500\}u=\-0\.06\\,\\mathcal\{R\}\_\{0\.29\}u\+0\.16\\,\\mathcal\{R\}\_\{1\.44\}u†\\dagger*Frac\. Burgers*— true:∂tu=0\.25​Dx1\.70​u−1\.00​u​ux\\partial\_\{t\}u=0\.25\\,D\_\{x\}^\{1\.70\}u\-1\.00\\,u\\,u\_\{x\}Frac\. Burgers0%∂tu=0\.25​Dx1\.70​u−1\.00​u​Dx1\.00​u\\partial\_\{t\}u=0\.25\\,D\_\{x\}^\{1\.70\}u\-1\.00\\,u\\,D\_\{x\}^\{1\.00\}u10%∂tu=0\.25​Dx1\.70​u−1\.00​u​Dx1\.00​u\\partial\_\{t\}u=0\.25\\,D\_\{x\}^\{1\.70\}u\-1\.00\\,u\\,D\_\{x\}^\{1\.00\}u
### 16Two\-dimensional settings and sensitivity

The two\-dimensional results use an example extension distributed in Online Resource 1\. It reuses the one\-dimensional implementations of the Gaussian test basis and the L1 Caputo adjoint, while its data generator includes a complex\-capable evaluator ofEα,1E\_\{\\alpha,1\}because directional Fourier multipliers are complex\. The generator and discovery code are independent in time: data use semi\-analytic Mittag–Leffler propagation, whereas discovery evaluates the Caputo target through the transposed L1 matrix\. The field contains five conjugate Fourier\-mode pairs, and the supplied archive records the data hashes, software environment, per\-seed estimates, validation curves, and complete summary\.

The default spatial width applies the one\-dimensional paper rule independently toxxandyy\. Table[22](https://arxiv.org/html/2608.12879#S16.T22)compares three fixed fractions of the domain length on Benchmark \([28](https://arxiv.org/html/2608.12879#S5.E28)\) at5%5\\%noise\. The inherited rule and0\.10​L0\.10Lrecover the complete structure in every seed, whereas0\.16​L0\.16Land0\.24​L0\.24Lselect only two terms\. This complements Appendix[12](https://arxiv.org/html/2608.12879#S12), where narrower localised Gaussian test\-function windows eventually degrade operator recovery: no two\-dimensional tuning was needed here, but test\-window scale remains problem dependent\.

Table 22:Spatial\-window sensitivity for Benchmark \([28](https://arxiv.org/html/2608.12879#S5.E28)\) at5%5\\%multiplicative noise \(five seeds\)\. Errors are conditioned on support and direction recovery\.Width ruleSupport/directionOperatoreαe\_\{\\alpha\}eβmaxe\_\{\\beta\}^\{\\max\}eξmaxe\_\{\\xi\}^\{\\max\}Inherited paper rule5/55/55/55/50\.00177±0\.000410\.00177\\pm 0\.000410\.00483±0\.000450\.00483\\pm 0\.000450\.01336±0\.001970\.01336\\pm 0\.001970\.10​L0\.10L5/55/55/55/50\.00150±0\.000340\.00150\\pm 0\.000340\.00875±0\.000800\.00875\\pm 0\.000800\.00669±0\.001190\.00669\\pm 0\.001190\.16​L0\.16L0/50/50/50/5–––0\.24​L0\.24L0/50/50/50/5–––Table[23](https://arxiv.org/html/2608.12879#S16.T23)compares Benchmark \([27](https://arxiv.org/html/2608.12879#S5.E27)\) at the reported and refined spatial grids\. Both resolutions recover all five noise levels and seeds\. The larger tensor\-product set of weak rows improves the coefficient error at20%20\\%noise while leaving the order errors of comparable scale; the result serves as a grid\-stability check\. Sparse observation patterns would require a different row\-construction and validation analysis and are outside the present scope\.

Table 23:Spatial\-resolution check for Benchmark \([27](https://arxiv.org/html/2608.12879#S5.E27)\)\. Recovery counts aggregate 25 runs; errors are at20%20\\%noise\.
### 17Comparison with a contemporary neural fractional\-discovery framework

The Strong\-Pareto versus Weak\-Pareto \(no polishing\) ablation isolates the candidate\-library effect under a common selector\. We also compare Weak\-Pareto with Yu et al\.[23](https://arxiv.org/html/2608.12879#bib.bib19)on the advection–diffusion benchmark\. The full Yu et al\. framework combines neural field reconstruction, automatic differentiation of integer derivatives, pointwise Gauss–Jacobi fractional derivatives, sparse regression, and global optimisation\. An optimiser\-only variant replaces the neural reconstruction with a deterministic quintic spline to isolate the downstream derivative and selection stages\.

The comparison shares the FADE field, nominal equation, target orders, noise realisations, seeds, recovery tolerances, and CPU environment\. It is not operator\-identical: Weak\-Pareto uses the periodic directional spectral operator that generated the field, whereas the Yu adaptation retains the one\-sided finite\-terminal Gauss–Jacobi approximation; the Riesz reaction–diffusion cases are outside its declared operator scope\. The adapter fixes the training–validation split, estimates coefficients and the STRidge penalty from training rows only, and uses deterministic seeds\. Its changes, budgets, and provenance controls are documented in Online Resource 1\. Because the upstream snapshot is not redistributed, byte identity remains unverified\. Runtime covers each complete fitting framework\. Table[24](https://arxiv.org/html/2608.12879#S17.T24)reports this controlled comparison\.

Table 24:Weak\-Pareto versus the adapted neural fractional\-discovery framework of Yu et al\.[23](https://arxiv.org/html/2608.12879#bib.bib19)on the advection–diffusion benchmark \(five seeds\)\. All rows use the same data, noise realisations, recovery tolerances, scoring convention, and CPU environment\. Errors are mean±\\pmstandard deviation over runs with correct support and power; runtime is mean±\\pmstandard deviation over all five seeds\. Table[12](https://arxiv.org/html/2608.12879#S5.T12)reports a separate uncached single\-seed timing\.Weak\-Pareto recovers the complete FADE operator in all five seeds at0%0\\%,1%1\\%, and5%5\\%noise, with order and coefficient errors of a few percent\. The adapted neural fractional\-discovery framework recovers1/51/5,2/52/5, and4/54/5operators; the non\-monotone counts reflect run\-to\-run variability, not evidence that noise improves recovery\. The optimiser\-only variant never recovers the complete operator: at0%0\\%its spatial\-order error is about0\.200\.20, and at positive noise it drops the fractional\-diffusion term\. On that CPU, Weak\-Pareto takes6\.66\.6–6\.76\.7s per run, versus279279–332332s for the adapted framework and1111–2525s for the optimiser\-only variant\. These are implementation\-level timings for the tested configurations\.

### Acknowledgements

The authors would like to thank Velmurugan Gandhi for his helpful discussions on fractional differential equations\.

### Supplementary material

Online Resource 1\.Source code, benchmark datasets, archived reference outputs, tutorials, tests, and scripts for reproducing the numerical results and figures reported in this article\. We will also maintain the code and reproducibility materials at[https://github\.com/Pongpisit\-Thanasutives/Weak\-Pareto](https://github.com/Pongpisit-Thanasutives/Weak-Pareto)\.

## References

- S\. L\. Brunton, J\. L\. Proctor, and J\. N\. KutzDiscovering governing equations from data by sparse identification of nonlinear dynamical systems\.Proceedings of the National Academy of Sciences113\(15\),pp\.3932–3937\.External Links:[Document](https://dx.doi.org/10.1073/pnas.1517384113)Cited by:[§1](https://arxiv.org/html/2608.12879#S1.p1.1)\.
- Cranmeret al\.\(2020\)M\. Cranmer, A\. Sanchez Gonzalez, P\. Battaglia, R\. Xu, K\. Cranmer, D\. Spergel, and S\. HoDiscovering symbolic models from deep learning with inductive biases\.InAdvances in Neural Information Processing Systems,Vol\.33,pp\.17429–17442\.Cited by:[§1](https://arxiv.org/html/2608.12879#S1.p5.1)\.
- Faselet al\.\(2022\)U\. Fasel, J\. N\. Kutz, B\. W\. Brunton, and S\. L\. BruntonEnsemble\-SINDy: robust sparse model discovery in the low\-data, high\-noise limit, with active learning and control\.Proceedings of the Royal Society A478\(2260\),pp\.20210904\.External Links:[Document](https://dx.doi.org/10.1098/rspa.2021.0904)Cited by:[§1](https://arxiv.org/html/2608.12879#S1.p4.1)\.
- Gulianet al\.\(2019\)M\. Gulian, M\. Raissi, P\. Perdikaris, and G\. KarniadakisMachine learning of space\-fractional differential equations\.SIAM Journal on Scientific Computing41\(4\),pp\.A2485–A2509\.External Links:[Document](https://dx.doi.org/10.1137/18M1204991)Cited by:[Table 1](https://arxiv.org/html/2608.12879#S1.T1.4.4.1.1),[§1](https://arxiv.org/html/2608.12879#S1.p5.1)\.
- Gurevichet al\.\(2019\)D\. R\. Gurevich, P\. A\. K\. Reinbold, and R\. O\. GrigorievRobust and optimal sparse regression for nonlinear PDE models\.Chaos: An Interdisciplinary Journal of Nonlinear Science29\(10\),pp\.103113\.External Links:[Document](https://dx.doi.org/10.1063/1.5120861)Cited by:[§1](https://arxiv.org/html/2608.12879#S1.p4.1)\.
- Kilbaset al\.\(2006\)A\. A\. Kilbas, H\. M\. Srivastava, and J\. J\. TrujilloTheory and applications of fractional differential equations\.Elsevier,Amsterdam\.Cited by:[§1](https://arxiv.org/html/2608.12879#S1.p2.1),[§2\.1](https://arxiv.org/html/2608.12879#S2.SS1.p1.2)\.
- Manganet al\.\(2017\)N\. M\. Mangan, J\. N\. Kutz, S\. L\. Brunton, and J\. L\. ProctorModel selection for dynamical systems via sparse regression and information criteria\.Proceedings of the Royal Society A473\(2204\),pp\.20170009\.External Links:[Document](https://dx.doi.org/10.1098/rspa.2017.0009)Cited by:[§1](https://arxiv.org/html/2608.12879#S1.p1.1)\.
- Messenger and Bortz \(2021a\)D\. A\. Messenger and D\. M\. BortzWeak SINDy for partial differential equations\.Journal of Computational Physics443,pp\.110525\.External Links:[Document](https://dx.doi.org/10.1016/j.jcp.2021.110525)Cited by:[Table 1](https://arxiv.org/html/2608.12879#S1.T1.4.2.1.1),[§1](https://arxiv.org/html/2608.12879#S1.p4.1)\.
- Messenger and Bortz \(2021b\)D\. A\. Messenger and D\. M\. BortzWeak SINDy: Galerkin\-based data\-driven model selection\.Multiscale Modeling & Simulation19\(3\),pp\.1474–1497\.External Links:[Document](https://dx.doi.org/10.1137/20M1343166)Cited by:[§1](https://arxiv.org/html/2608.12879#S1.p4.1)\.
- Metzler and Klafter \(2000\)R\. Metzler and J\. KlafterThe random walk’s guide to anomalous diffusion: a fractional dynamics approach\.Physics Reports339\(1\),pp\.1–77\.External Links:[Document](https://dx.doi.org/10.1016/S0370-1573%2800%2900070-3)Cited by:[§1](https://arxiv.org/html/2608.12879#S1.p2.1)\.
- Panget al\.\(2019\)G\. Pang, L\. Lu, and G\. E\. KarniadakisfPINNs: fractional physics\-informed neural networks\.SIAM Journal on Scientific Computing41\(4\),pp\.A2603–A2626\.External Links:[Document](https://dx.doi.org/10.1137/18M1229845)Cited by:[§1](https://arxiv.org/html/2608.12879#S1.p5.1)\.
- Podlubny \(1999\)I\. PodlubnyFractional differential equations\.Academic Press,San Diego\.Cited by:[§1](https://arxiv.org/html/2608.12879#S1.p2.1),[§2\.1](https://arxiv.org/html/2608.12879#S2.SS1.p1.2)\.
- Reinboldet al\.\(2020\)P\. A\. K\. Reinbold, D\. R\. Gurevich, and R\. O\. GrigorievUsing noisy or incomplete data to discover models of spatiotemporal dynamics\.Physical Review E101\(1\),pp\.010203\.External Links:[Document](https://dx.doi.org/10.1103/PhysRevE.101.010203)Cited by:[§1](https://arxiv.org/html/2608.12879#S1.p4.1)\.
- Rudyet al\.\(2017\)S\. H\. Rudy, S\. L\. Brunton, J\. L\. Proctor, and J\. N\. KutzData\-driven discovery of partial differential equations\.Science Advances3\(4\),pp\.e1602614\.External Links:[Document](https://dx.doi.org/10.1126/sciadv.1602614)Cited by:[§1](https://arxiv.org/html/2608.12879#S1.p1.1)\.
- Satopaaet al\.\(2011\)V\. Satopaa, J\. Albrecht, D\. Irwin, and B\. RaghavanFinding a “kneedle” in a haystack: detecting knee points in system behavior\.In2011 31st International Conference on Distributed Computing Systems Workshops,Minneapolis, MN, USA,pp\.166–171\.External Links:[Document](https://dx.doi.org/10.1109/ICDCSW.2011.20)Cited by:[§4\.3](https://arxiv.org/html/2608.12879#S4.SS3.p2.1)\.
- Schaeffer and McCalla \(2017\)H\. Schaeffer and S\. G\. McCallaSparse model selection via integral terms\.Physical Review E96\(2\),pp\.023302\.External Links:[Document](https://dx.doi.org/10.1103/PhysRevE.96.023302)Cited by:[§1](https://arxiv.org/html/2608.12879#S1.p4.1)\.
- Schaeffer \(2017\)H\. SchaefferLearning partial differential equations via data discovery and sparse optimization\.Proceedings of the Royal Society A473\(2197\),pp\.20160446\.External Links:[Document](https://dx.doi.org/10.1098/rspa.2016.0446)Cited by:[§1](https://arxiv.org/html/2608.12879#S1.p1.1)\.
- Schmidt and Lipson \(2009\)M\. Schmidt and H\. LipsonDistilling free\-form natural laws from experimental data\.Science324\(5923\),pp\.81–85\.External Links:[Document](https://dx.doi.org/10.1126/science.1165893),[Link](https://www.science.org/doi/abs/10.1126/science.1165893)Cited by:[§1](https://arxiv.org/html/2608.12879#S1.p5.1)\.
- Stephany and Earls \(2024\)R\. Stephany and C\. EarlsWeak\-PDE\-LEARN: a weak form based approach to discovering PDEs from noisy, limited data\.Journal of Computational Physics506,pp\.112950\.External Links:[Document](https://dx.doi.org/10.1016/j.jcp.2024.112950)Cited by:[Table 1](https://arxiv.org/html/2608.12879#S1.T1.4.3.1.1),[§1](https://arxiv.org/html/2608.12879#S1.p4.1)\.
- Storn and Price \(1997\)R\. Storn and K\. PriceDifferential evolution – a simple and efficient heuristic for global optimization over continuous spaces\.Journal of Global Optimization11\(4\),pp\.341–359\.External Links:[Document](https://dx.doi.org/10.1023/A%3A1008202821328)Cited by:[§4\.2](https://arxiv.org/html/2608.12879#S4.SS2.p1.5)\.
- Tanget al\.\(2023\)M\. Tang, W\. Liao, R\. Kuske, and S\. H\. KangWeakIdent: weak formulation for identifying differential equation using narrow\-fit and trimming\.Journal of Computational Physics483,pp\.112069\.External Links:[Document](https://dx.doi.org/10.1016/j.jcp.2023.112069)Cited by:[Table 1](https://arxiv.org/html/2608.12879#S1.T1.4.2.1.1),[§1](https://arxiv.org/html/2608.12879#S1.p4.1)\.
- Thanasutiveset al\.\(2024\)P\. Thanasutives, T\. Morita, M\. Numao, and K\. FukuiAdaptive uncertainty\-penalized model selection for data\-driven PDE discovery\.IEEE Access12,pp\.13165–13182\.External Links:[Document](https://dx.doi.org/10.1109/ACCESS.2024.3354819)Cited by:[§1](https://arxiv.org/html/2608.12879#S1.p4.1)\.
- Yuet al\.\(2025\)X\. Yu, H\. Xu, Z\. Mao, H\. Sun, Y\. Zhang, Z\. Chen, D\. Zhang, and Y\. ChenA data\-driven framework for discovering fractional differential equations in complex systems\.Nonlinear Dynamics113\(18\),pp\.24557–24577\.External Links:[Document](https://dx.doi.org/10.1007/s11071-025-11373-z)Cited by:[Table 1](https://arxiv.org/html/2608.12879#S1.T1.4.5.1.1),[§1](https://arxiv.org/html/2608.12879#S1.p5.1),[Table 24](https://arxiv.org/html/2608.12879#S17.T24),[§17](https://arxiv.org/html/2608.12879#S17.p1.1),[§5\.4](https://arxiv.org/html/2608.12879#S5.SS4.p3.1),[§5\.7](https://arxiv.org/html/2608.12879#S5.SS7.p1.1)\.

### Statements and Declarations

Funding\.Pongpisit Thanasutives is supported by the research fund from the Special Postdoctoral Researcher \(SPDR\) program at RIKEN, Japan\.

Competing interests\.The authors have no relevant financial or non\-financial interests to disclose\.

Author contributions\.Pongpisit Thanasutives designed the method and software, performed the investigation and validation, analysed the results, and wrote the manuscript\. Yoshinobu Kawahara supervised the work and revised the manuscript\. Both authors reviewed and approved the final manuscript\.

Use of generative AI\.During manuscript preparation, the authors used OpenAI ChatGPT and Anthropic Claude for language revision, manuscript drafting, and code review\. The authors reviewed and verified all mathematical arguments, software implementations, results, and final text and take full responsibility for the work\.

Data and code availability\.Data and source code are provided in the Supplementary Material and will be maintained at[https://github\.com/Pongpisit\-Thanasutives/Weak\-Pareto](https://github.com/Pongpisit-Thanasutives/Weak-Pareto)\. Third\-party frozen\-soil data are not redistributed; access instructions and provenance are documented in the cited paper\.