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

The paper Robust data-driven discovery of fractional differential equations via weak formulations and Pareto-based subset selection by Pongpisit Thanasutives and Yoshinobu Kawahara proposes Weak-Pareto, a data-driven discovery framework for fractional partial differential equations (FPDEs) that combines an adjoint-consistent weak formulation of fractional terms with Pareto-based subset selection over discrete term types and continuous fractional orders.

The authors note that fractional equation discovery introduces two difficulties beyond the integer-order setting. First, fractional differentiation can amplify high-frequency measurement noise—for periodic spectral operators, an operator of order β scales a Fourier mode of wavenumber κ by a factor whose magnitude grows as κβ. Second, the unknown orders are continuous—a fixed dictionary must either use a coarse order grid (creating discretisation bias) or a dense grid (producing many nearly collinear columns that destabilise support selection).

Existing fractional-discovery methods still evaluate fractional derivatives pointwise, either directly on a mesh or after reconstructing the field with a neural network. The authors state that 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.

The paper claims three core contributions:

1. An adjoint-consistent weak library for fractional operators: "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. For nonlinear terms, the weak row is an averaged projection of the corresponding strong feature."

2. A branch-aware, continuous-order Pareto search: "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."

3. Controlled analyses and baseline comparisons: Including library-only ablations, fixed-dictionary ablations, and a controlled comparison assesses an adapted neural fractional-discovery framework.

Proposition 1 shows that for independent, zero-mean noise with variance σ2, the weak-feature perturbation satisfies Var⟨η, ω⟩ h = O(h t h x), meaning its variance is O((n t n x)−1) under grid refinement. In contrast, for a periodic pointwise feature, Var(X β η) ij ∼ (π(2β)/(2β+1))σ2h x(−2β), meaning "the variance of every unregularised positive-order strong-form feature (β > 0) diverges under spatial refinement, at a rate that increases with β."

Corollary 2 extends the variance result to multiplicative measurement noise, showing the conditional perturbation variance is O(h t h x).

Remark 1 derives the bias of nonlinear weak features: for p = 1 under additive i.i.d. noise, the discrete-adjoint feature has bias σ2h t h x tr(A* h D φk). For the periodic directional multiplier of order β, the leading bias behaves as σ2cos(πβ/2)(πβ/(β+1))h x(−β)∫ Q φ k.

Remark 2 analyzes noise in the Caputo target, noting that the same noisy observation u(0, x j) enters every temporal contribution through the Caputo initial-value subtraction in the subunit branch, and that in the superunit branch, the transpose of the composed L1–first-difference operator places a large opposite-sign pair of weights on the first two time samples.

Proposition 3 characterizes local identifiability of spatial fractional orders: after coefficient refitting, the smallest first-order change in the model mean is ξ j2‖(I − P Θ)θ̇ j‖2. This motivates the diagnostic S βj = ‖(I − P Θ)θ̇ j‖2/‖θ j‖2.

The model class is the parsimonious fractional equation T m α,αu = Σ ξ j u(p j) X β ju. Each candidate is represented by the tuple M = (m α, α, p, β, ξ), where m α records the temporal branch (subunit for 0 < α < 1, exact integer at α = 1, superunit for 1 < α < 2).

The weak formulation uses test functions φ k(t,x) = ϑ l t(t)ψ l x(x), producing scalar equations ⟨T m α,αu, φ k⟩ = Σ ξ j⟨u(p j)X β ju, φ k⟩. The construction rests on the adjoint identity ⟨Lf, φ k⟩ = ⟨f, L*φ k⟩ + B L(f, φ k).

For the Caputo target, fractional integration by parts gives ⟨C 0D t α u, φ k⟩ = ⟨u − P n−1,0u, (tD T α)φ k⟩. For linear spatial terms, θ k,β = ⟨X βu, φ k⟩ = ⟨u, X βφ k⟩. For nonlinear terms, ⟨u pX βu, φ k⟩ = ⟨u, X β(u pφ k)⟩.

The numerical library uses the discrete adjoint A* h = W−1A TW, satisfying ⟨Af, φ⟩ h = ⟨f, A* hφ⟩ h.

The search procedure solves best-subset problems indexed by support size c, using differential evolution for continuous orders and ridge regression for coefficients. The elbow selection uses a signed-chord criterion on the benefit–complexity front, where the selected interior point maximises its signed vertical excess above the chord joining the first and last points.

Four periodic main benchmarks are used: FADE (D t 0.8u = −∂ xu + 0.5D x 1.7u), two Riesz reaction–diffusion equations, and a nonlinear fractional Burgers equation (∂ tu = −u∂ xu + 0.25D x 1.7u).

At 10% multiplicative noise, 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 in 3/5 space-fractional runs and all five time–space runs, but no run satisfies the full operator criterion at 10% noise.

At 0% noise, both Riesz benchmarks recover the complete operator in all five seeds. The time–space case reaches 2/5 operator recoveries at 2% noise; no positive-noise row for the space-fractional case does.

For the semi-analytic equation C 0D t 1.65u = 0.12D x2u, Both methods recover the clean superunit operator in all five runs. At 0.5% noise, Weak-Pareto retains 5/5 branch and operator recovery, while the strong-form comparator selects the wrong temporal branch in every run. At 1%, Weak-Pareto still selects the superunit branch in all five runs, but four estimates reach the upper search bound α = 1.85, leaving only 1/5 complete operator recoveries.

Under multiplicative noise, Weak-Pareto recovers the correct support in all five seeds at every tested level up to 20% on both benchmarks (FADE and fractional Burgers). The strong-form framework does not recover the correct Burgers support in any noisy run and in all but the 5% FADE condition. Under additive Gaussian noise, Weak-Pareto recovers the correct support in all five runs for both equations and the complete operator in 4/5 seeds for each; the strong-form framework does not recover the correct support in any run for either equation.

The library-only contrast is decisive: With the same continuous-order selector and no polishing, Strong-Pareto achieves no complete FADE operator recovery at 10% noise, whereas Weak-Pareto recovers all five. The fixed-dictionary baseline (Weak Grid-STRidge) does not recover the correct FADE support in any of the five runs even though its fixed grid contains nodes within 0.01 of both true orders. Disabling exact-order polishing leaves recovery unchanged and changes e β max only from 0.08 to 0.07.

The encoding is augmented with a direction label d j ∈ x, y. Both benchmarks retain complete direction and operator recovery throughout the noise sweep across 25 runs per benchmark (five seeds × five noise levels). At 20% noise, errors remain small: for Benchmark (27), (e α, e β max, e ξ max) = (0.00081 ± 0.00054, 0.00217 ± 0.00140, 0.0111 ± 0.0023).

On the advection–diffusion benchmark, Weak-Pareto recovers the complete FADE operator in all five seeds at 0%, 1%, and 5% noise, with order and coefficient errors of a few percent. The adapted neural fractional-discovery framework recovers 1/5, 2/5, and 4/5 operators. Runtime: Weak-Pareto takes 6.6–6.7 s per run, versus 279–332 s for the adapted framework.

For the fractional Kelvin model applied to irregular experimental data, The silt estimate reproduces the reference order and time scale closely. For clay, the order differs by 0.098 and t ret by about 78%, indicating substantial parameter uncertainty despite a small integral residual.

The authors identify two empirical limits: In the Riesz reaction–diffusion experiments, support and powers are often retained while nearby spatial orders and small reaction coefficients remain difficult to distinguish. And in the fixed-support superunit diagnostic, the correct branch is retained at 1% noise, but four of five temporal-order estimates reach the upper search bound.

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.

Improvement 1: Noise-Robust PDE Discovery via Adjoint-Consistent Weak Formulation

The AI system can now discover fractional PDEs from noisy spatiotemporal data without amplifying high-frequency noise. By integrating the adjoint-consistent weak library, it avoids pointwise differentiation of measured fields for linear terms, reducing variance from O(h x(-2β)) (diverging under refinement) to O(h t h x) (vanishing under refinement). The system can recover operators (e.g., FADE, fractional Burgers) at up to 20% multiplicative noise across all seeds, whereas strong-form baselines fail entirely. It can also handle multiplicative noise with conditional perturbation variance O(h t h x), enabling robust discovery in experimental settings.

Improvement 2: Continuous-Order Pareto Search with Branch-Aware Encoding

The system can jointly optimize discrete term types (integer powers) and continuous fractional orders without discretization bias. It uses differential evolution over orders and analytic ridge regression for coefficients, treating subunit (0<α<1), integer (α=1), and superunit (1<α<2) Caputo branches as distinct. This eliminates the need for dense order grids that cause collinearity, and it correctly identifies superunit temporal orders (e.g., α=1.65) even at 1% noise, where strong-form methods select the wrong branch in every run.

Improvement 3: Identifiability Diagnostics for Fractional Orders

The system can quantify local identifiability of each spatial fractional order using the diagnostic S βj = (I − P Θ)θ̇ j2/θ j2. This allows it to flag orders that are poorly constrained by data (e.g., in Riesz reaction–diffusion cases where nearby orders are indistinguishable). The AI can then either report uncertainty or adaptively refine the search around identifiable orders, improving reliability in underdetermined scenarios.

Improvement 4: Bias-Aware Nonlinear Feature Handling

For nonlinear terms (e.g., u p X β u), the system uses averaged projections of strong features, but it can now compute and report the leading bias (e.g., σ2h t h x tr(A* h D φk) for p=1, or σ2cos(πβ/2)(πβ/(β+1))h x(-β)∫φ k for periodic multipliers). This enables the AI to correct or down-weight biased nonlinear terms during discovery, improving accuracy in advection–diffusion and Burgers-type equations where nonlinearity dominates.

Improvement 5: Efficient Neural-Baseline Replacement

The system can replace neural-network-based fractional PDE discovery with a weak-formulation approach that runs 40–50× faster (6.6–6.7 s vs. 279–332 s per run) while achieving higher operator recovery rates (5/5 vs. 1/5–4/5 at various noise levels). This makes real-time or large-scale discovery feasible on standard hardware, and it can be deployed in iterative experimental design loops.

Improvement 6: 2D Direction-Aware Discovery

By augmenting the encoding with direction labels (d j ∈ x, y), the system can discover anisotropic fractional PDEs in two spatial dimensions, recovering both directions and operators across all noise levels (0–20%) with errors below 0.003 in orders and 0.012 in coefficients. This extends applicability to realistic transport and reaction–diffusion systems.

Improved AI System Capabilities

• Discover fractional PDEs from noisy, irregular, or sparse data (e.g., frozen-soil creep experiments) with quantified parameter uncertainty.

• Automatically select between temporal branches (subunit/integer/superunit) and spatial operators without user priors.

• Provide confidence intervals for fractional orders and coefficients, flagging non-identifiable parameters.

• Handle both additive and multiplicative noise, with theoretical guarantees on variance reduction.

• Scale to 2D problems and run in seconds, enabling interactive model discovery in laboratory or field settings.

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 contemporary neural baseline.