Constrained Deep Operator Networks with a Fully Discrete Error Bound for a Material-Coordinate Fractional-Diffusion Benchmark

In the original timing run, batched CP-DON inference took 11.53 μs per case against 711.47 μs per 32×20 L1 solve, a ratio of 61.7. This is an inference-only benchmark under the stated CPU setting. Table 6 shows that the neural inference cost of DON and C-DON was almost the same. The comparison does not establish the speed-up of a certified CP-DON workflow, because σmin(A) computation and residual-certificate assembly were excluded. It also compares a batched forward pass over 80 cases with a sequential per-case Python solver, so the two sides do not share the same batching, threading or input-preparation assumptions. A re-measurement on documented hardware (Table 6, lower block) gave 16.7 ± 5.9 μs per case for batched CP-DON and 1270 μs per solve, a ratio of about 76. For a single-case call (batch size 1), however, the network took about 800 μs, only 1.6 times faster than the solver. Computing the certificate with an explicit SVD of the 570×570 system took about 68 ms per case, roughly 50 times the cost of a solve. On this benchmark the certified workflow is therefore slower than solving directly. The timing claim is limited to amortised, batched neural inference on the tested grid.

3.5. Training-Seed Robustness and Equal-Epoch Comparison

To address the single-seed limitation, all three models were retrained with five additional seeds (1–5; different initialisation and mini-batch order, identical datasets) using the settings of Table 2 (140 epochs for DON and C-DON, 170 for CP-DON). CP-DON was additionally retrained with the 140-epoch cap used for the other two models. Early stopping did not trigger in any retraining run. Table 7 reports the mean ± SD over seeds of the test-set mean relative L2 error.

The ordering DON > C-DON > CP-DON (170 epochs) held in every seed and every regime. The reductions of CP-DON relative to DON (61.2 ± 4.3%, 56.3 ± 3.6% and 41.5 ± 4.0%) agree with the single-seed values of 58.4%, 57.9% and 41.7%, which lie within one standard deviation. The trial map accounted for 82 ± 9%, 80 ± 7% and 82 ± 4% of the CP-DON improvement over DON in the three regimes, consistent with the single-seed 87%.

The additional benefit of the residual term is less secure. At the 170-epoch budget CP-DON lowered the error relative to C-DON by 21.2 ± 6.0%, 19.8 ± 4.3% and 11.6 ± 2.6%. When limited to the same 140 epochs, however, CP-DON was better than C-DON in only 2 of 5 seeds in every regime (mean relative difference +5.4%, +4.0% and +1.8%, i.e., slightly worse on average). The gain attributed to the residual term therefore cannot be separated from the longer and costlier training; on the hardware a CP-DON run took about 55–70 s, against about 10 s for C-DON and 9 s for DON (seeds 3–5, uncontended). An equal-wall-clock comparison was not performed. Finally, the seed-averaged final-time geometry errors were 0.053 (DON), 0.060 (C-DON) and 0.056 (CP-DON), so the constraints did not improve the final-time error under geometry shift.

3.6. Discretisation Convergence and Continuum Context

Convergence of the L1 scheme was examined on nested grids with (Nx−1, Nt−1) = (31k, 19k) for k = 1, 2, 4, 8, 16, 32, so that every grid contains the 32×20 benchmark nodes (k = 1 reproduced the dataset solutions to within 6×10⁻⁸). For the 60 audited parameter cases the observed order, estimated at the common nodes from the differences between three successive grid levels (k, 2k, 4k), increased from 0.93 at k = 1 to 0.99 at k = 8 (case range 0.81–1.07), i.e., first-order convergence under joint refinement (Table 8). A constant-order, constant-geometry special case with initial profile sin(πy) has the exact solution u(y,t) = sin(πy)Eα(−κπ²L⁻²tα), where Eα is the Mittag-Leffler function and tα denotes t raised to the power α. For α = 0.4, 0.6, 0.8 and 0.95 the same scheme converged to it with observed orders between 0.94 and 1.08, and the relative L2 error fell from 1.4–2.7% on the 32×20 grid to 0.05–0.08% on the finest grid. The scheme therefore converges to the continuum solution at first order, but the 32×20 benchmark grid is coarse.

Relative to the finer reference, the error of CP-DON is 1.8–2.8 times larger than relative to the 32×20 discrete reference (Table 8), because the training data themselves carry the discretisation error. The ranking of the three models is unchanged, so the comparative conclusions stand; what changes is the meaning of the absolute numbers. They measure how well the surrogate reproduces the 32×20 solver, not the continuum PDE. In the audited cases the discretisation error δ, measured in the same ℓ2 norm as the certificate, was 6–38% of the certificate (mean 16%) and exceeded the surrogate error in 49 of 60 cases. A combined indicator η + δ bounds the surrogate error against the finer reference by the triangle inequality, but δ is itself estimated from a finer numerical solution and is not a rigorous quantity.

4. DISCUSSION

The results support a practical but limited claim. For this material-coordinate benchmark, a constrained operator surrogate, a discrete residual loss and a residual/singular-value bound can be combined in one pipeline. The evidence is confined to the selected discrete model and parameter family, to a training procedure whose variability was quantified with only five additional seeds (Section 3.5), and to a discretisation whose own error (1.7–2.0%) exceeds the surrogate error measured against it (Section 3.6).

4.1. Exact Admissibility as an Effective Operator Bias

The main conclusion is that structural encoding was not ‘cosmetic’. C-DON lowered the in-distribution error by 51.1% relative to DON without any physics residual. In absolute terms, the trial map and the residual term contributed 0.0088 and 0.0013, respectively, about 87% and 13% of the 0.0101 DON-to-CP-DON improvement in the single-seed run. Across five further seeds the trial-map share was 82 ± 9% (80 ± 7% and 82 ± 4% under the two shifts), so the dominance of the structural component is robust, whereas the 13% residual share is uncertain (Section 4.2). The constrained models also showed no initial, boundary or sign violations in the measured data, except for the endpoint caveat in Section 2.3. The result is consistent with exact-condition approaches [26], which restrict the output class and can ease learning because the network need not approximate properties already known for the benchmark.

The result also warns against rating a scientific surrogate on trajectory error alone. The average error for DON was low, but the errors were negative at about 3% of the grid points and non-zero on the boundaries. These violations may be intolerable if u is a concentration or a probability-like value. A soft-penalty alternative would be to assign weights to initial and endpoint losses and positivity losses, and would possibly leave behind some violations. The exact map maps these properties from optimisation objectives to identities.

This restriction is fundamental, not incidental. The trial map (6) is too restrictive when the true solution changes sign, when boundaries are non-homogeneous or uncertain, when an interior zero is later filled by diffusion, or when the solution arises from a zero initial condition with a source term. The map should therefore not be applied blindly; a trial representation should instead be built to match the known admissible class.

4.2. Incremental Contribution of Physics Regularisation

The residual term added a smaller benefit once the exact constraints were imposed. In the single-seed run CP-DON lowered the error relative to C-DON by 15.0% in distribution, 16.5% under fractional-order shift and 6.7% under geometry shift, and it also lowered the normalised residual in all regimes. Over five further seeds, at the 170-epoch budget, the reductions were 21.2 ± 6.0%, 19.8 ± 4.3% and 11.6 ± 2.6%. This is not an equal-budget result: when CP-DON was limited to the 140 epochs given to C-DON, it was better than C-DON in only 2 of 5 seeds in every regime (Table 7), and a CP-DON epoch was about six times as expensive as a C-DON epoch. The residual benefit is therefore not established independently of the longer, costlier training. The trend resembles the PINO and VINO studies [22, 23], where equation information improves data efficiency or high-resolution fidelity, but the gain is modest here because the training labels are exact and λphys is small.

The ablation clarifies a point that is often confused: a physics loss and hard constraints address different problems. The trial map guarantees admissibility, whereas the residual term promotes consistency with the discrete evolution law. The normalised residual of C-DON was slightly larger than that of DON in distribution and under fractional-order shift (0.280 versus 0.256 and 0.310 versus 0.279), even though C-DON is admissible and more accurate; a small residual therefore does not imply admissibility, and admissibility does not imply a small residual. Where both are affordable, combining them is reasonable, but the present evidence supports the combination only at the cost of the longer training noted above.

The use of the discrete L1 residual also has consequences. It avoids automatic differentiation of a fractional derivative, which is costly and noisy, and it gives a training objective that matches the data generator. The regulariser is therefore a residual regulariser tied to the chosen 32×20 L1 discretisation, not to the continuum fractional PDE. The resulting CP-DON learns a surrogate aligned with this discrete operator and not an abstract continuum operator. Resolution transfer was not tested, and the model is not discretisation invariant.

4.3. A Rigorous but Discrete and Conservative Certificate

The reliability result is deterministic and discrete. All audited predictions satisfied the residual/singular-value bound. The certificate uses only the assembled matrix, the right-hand side of the discrete system and the predicted trajectory, and it needs no reference solution. It is valid for parameter values outside the training ranges provided that the same fully discrete linear system is correctly assembled and is nonsingular. The role is conceptually similar to residual-based a posteriori estimation for variational PINNs [29] and to variational neural approximation with error control [33], but the bound here is specifically a stability bound for the assembled discrete linear system and, unlike those estimators, controls no continuum quantity.

Three qualifications are required. First, the certificate bounds the error with respect to the discrete reference z*, not the continuum solution. The nested-grid study (Section 3.6) shows that the 32×20 solution differs from a 992×608 reference by 1.7–2.0% on average (up to 5.8%), and this discretisation error exceeded the surrogate error in 49 of 60 audited cases. A continuum certificate would need an additional discretisation estimate, for example by Richardson extrapolation, which the observed first-order convergence would support but which was not pursued. Second, with efficiency indices of 1.9–25.4 (case-wise medians 8.9–15.7), the σmin-based bound is conservative: a safety bound rather than an exact estimator. Third, the explicit computation of σmin(A) is expensive (about 68 ms per case, roughly 50 times the cost of a solve on this benchmark) and was excluded from the speed-up. At large dimension a low-cost certified lower bound, an energy-norm estimate, a problem-specific stability constant or a preconditioner would be required. The certificate provides no guarantee on the continuum PDE error, modelling error, unmodelled physics, inaccurate numerical assembly, floating-point perturbations, different topologies, or a scalable production certificate.

4.4. Extrapolation, Practical Significance, and Limitations

Fractional-order extrapolation was easier than geometry extrapolation. The relative L2 error was 1.00% for the order-shift set, and 1.92% for unseen interval lengths and rates, for CP-DON. The geometry error at the end of the run was greater than 5% and the amplitude bias is seen in Fig. (7) without any change in the qualitative mode. This behavior is consistent with the other literature on geometry-operator mapping across a domain change [18, 21]. While mapping across a changing domain is feasible, it is a different matter to extrapolate beyond the geometric support. The present tests are extrapolations of parameter ranges within a single PDE family, dimension, topology, affine domain map, discretisation and initial-profile family. They are not evidence of general out-of-distribution robustness, which would require changes of family rather than of parameter range.

In practical terms, the framework is most believable as an amortised proxy for multiple queries in a parameter range within a characterised range; the certificate as a screening device for the case level. The novelty of this reliability component is the incorporation of the discrete bound in this screening workflow, rather than the presentation of the underlying singular-value inequality as novel. The best deployment model is accept-or-solve: the surrogate predicts, the discrete certificate is evaluated, and the numerical solver is applied when the certificate is above a certain tolerance. This hybrid policy is more justifiable than the pure replacement of autonomous solvers. The reported 61.7-fold speed-up is an inference-only figure under the specified setting: it excludes training, data generation and certificate cost, it compares batched neural evaluation with a per-case solver, and it does not reflect the much larger systems of realistic multidimensional models. In the re-measurement the explicit SVD certificate cost far more than the solve itself (Section 3.4), so the accept-or-solve policy is practical only with a cheaper certificate.

Several limitations guide follow-up work. The study is a one-dimensional, linear and synthetic study. The order only changes with time; the motion is interval affine, and there is no topology change. The starting profiles in the first family make a smooth low dimensional family. Training variability was quantified with five additional seeds on fixed datasets; dataset resampling was not repeated, and the bootstrap intervals in Table 3 refer to test-case variability only. The physics weight and the longer CP-DON training horizon were chosen in a validation/pilot process; the equal-epoch retraining (Table 7) showed that the residual benefit depends on this longer horizon, and an equal-wall-clock comparison was not performed. The discretisation error of the benchmark (1.7–2.0%) is comparable to or larger than the surrogate error. Singular-value decompositions were expensive, so a systematic subset of 60 cases was used for the evaluation of the residual certificate. Last, but not least, no comparison was made with FNO, geometry-specific graph operators or adaptive numerical solvers, which were chosen for the study to isolate the proposed components and not to make them state-of-the-art.

Future research directions involve extending the trial-map concept to source-driven and sign-changing problems, to nonlinear maps and 2-D geometries, to report training variability across cases, to add case-wise certificate-efficiency dispersion, to replace explicit SVD with a scalable stability lower bound, and to add a discretisation-error estimate for a continuum-aware certificate. There should also be a multi-resolution training before the claim of continuum or mesh transfer.