Papers

Surrogate models deployed on real physical systems rarely see data at fixed resolutions: sensor configurations vary across deployments and may evolve over time as sensing infrastructure changes. We introduce the Latent Twin Operator (LTO), a latent-space surrogate for time-evolving PDEs with a Convolutional Conditional Neural Process (ConvCNP)-style encoder and decoder that accepts a context set of $N$ sensor observations with arbitrary placement and can be queried at any resolution. A learned latent evolution map advances the encoded state directly between arbitrary time points, decoupling temporal evolution from the observation and query discretizations. We derive an explicit $\mathcal{O}(N^{-2/(3D)})$ rate for the context discretization error in spatial dimension $D$ under quasi-uniform refinement. We verify the predicted decay empirically on a 2D heat-equation benchmark. Across time-dependent PDE benchmarks, LTO achieves strong accuracy under one-step comparisons and transfers from native to coarser spatial resolutions with fixed parameters. On Navier–Stokes, its direct latent evolution further reduces error over longer prediction horizons relative to recursive evaluation.

Many learning tasks map an input distribution to an output distribution. A natural way to model such an operator is to transform each input sample using a continuous function that may depend on the entire input distribution, and then take the distribution of the transformed samples. This defines a measure-dependent pushforward model and includes measure-theoretic formulations of transformers. We ask when such models can approximate arbitrary continuous operators between spaces of probability measures. We first show that universal approximation fails when atomic inputs are allowed: some continuous measure-to-measure operators that split or redistribute atomic mass cannot be approximated arbitrarily well by deterministic pushforward models. We then introduce the uniform level set condition, which requires a continuous measure-dependent scalarization whose shrinking level set neighborhoods carry uniformly vanishing mass over the input family. This condition is satisfied, in particular, by compact families of absolutely continuous measures. On every compact family satisfying this condition, we prove that any continuous measure-to-measure operator with outputs of finite $p$-th moment can be uniformly approximated, in the $p$-Wasserstein distance, by continuous measure-dependent pushforwards. Combining our theorem with existing approximation results for measure-dependent in-context maps yields universal approximation by measure-theoretic transformers. We also extend the framework to continuously-varying source measures, yielding a corresponding universality result for a class of pushforward models that are closely aligned with cross-attention architectures.

We construct the first convex planar counterexample to the Schiffer and Pompeiu conjectures: a bounded strictly convex non-disc domain with real-analytic boundary carrying a nonconstant solution of $\Delta u+u=0$ with $u=1$ and $\partial_\nu u=0$ on the boundary. Together with the higher-dimensional constructions, this gives convex non-ball counterexamples in every dimension from two to eighteen, and in dimensions twenty and twenty-one, with real-analytic boundaries diffeomorphic to spheres and indicator Fourier transforms vanishing on the unit sphere. In dimensions $4$, $6$, $8$, $10$ and $14$, planar reductions through compact Lie algebras of rank two use Harish-Chandra's radial part formula and Kostant's convexity theorem. In dimensions $3$, $5$ and $7$, an axisymmetric formulation on the unit ball gives a quartic equation in weighted coefficient spaces. The exact two-sided inverse of the residual-corrected linearisation combines a Dirichlet Helmholtz inverse, a holomorphic boundary real-part problem and the material-derivative identity; it gives the planar Schiffer domain and, by one dimension-parametrised argument, the domains in dimensions $5$, $9$, $11$ to $13$, $15$ to $18$, $20$ and $21$. We also construct the first convex planar counterexample to Berenstein's conjecture: a strictly convex non-disc domain with real-analytic boundary carrying a sign-changing solution of $\Delta u+u=0$ with $u=0$ and $\partial_\nu u=1$ on the boundary. Planar existence for this problem is due to Colbrook, Sadeghi and Stepaniants. All existence proofs in this paper are computer-assisted.

The evaluation cost of transformers with softmax attention scales quadratically with sequence length. Kernel attention addresses this by replacing softmax with a more general kernel function. In this paper, we aim to identify kernels that retain the expressivity of attention while enabling quasi linear computation. To quantify expressivity, we introduce a capacity for each kernel, measuring the maximum sequence length for which the attention matrix can approximate the identity. A higher capacity thus indicates greater expressivity. We show that expressive kernels like softmax, Gauss, and Laplace have infinite capacity. In contrast, common quasi linear kernels, such as those derived from finite dimensional feature maps, exhibit finite capacity. As a solution, we propose additive kernels constructed from univariate spline and polynomial exponential kernels. We prove that these maintain infinite capacity while allowing quasi linear computation via sorting. Finally, we implement additive sorting kernels efficiently and benchmark them against modern softmax backends, demonstrating advantages for long sequences.

For systems with steep gradients, sharp interfaces, or severe spatio-temporal coupling, Physics-informed neural networks (PINNs) suffer from spectral bias, geometric inflexibility, and boundary constraint conflicts, which undermine accuracy and convergence. To overcome these issues, we propose a geometry-adaptive and constraint-enhanced PINN (GAC-PINN). The framework comprises four components: a gradient-driven adaptive grid mapping (AGM) for diffeomorphic point concentration with Jacobian regularization, an adaptive bandwidth hard-constraint ansatz with spatially-varying boundary transition widths, a Gaussian Fourier feature mapping as a spectral preconditioner to further enhance high-wavenumber representation, and an operator-aware router that automatically selects the appropriate hard-constraint construction based on whether the governing PDE contains temporal derivatives. An AGM callback mechanism and a three-stage training strategy ensure stable coordination. Benchmarks including the viscous Burgers equation, a sharp-peaked 2D Poisson problem, and the Allen-Cahn phase-transition equation show that GAC-PINN attains relative (L^2) errors of ((1.747\pm 0.450)\times 10^{-4}), ((2.868\pm 0.947)\times 10^{-5}), and ((1.756 \pm 0.712)\times 10^{-3}), respectively, consistently outperforming the baselines. Ablation studies further reveal that AGM alone yields a substantially lower error than residual-based adaptive refinement (RAR), while RAR becomes beneficial only when combined with FFM, demonstrating a context-dependent module interaction. Convergence analysis verifies rapid error reduction and saturation with increasing resolution, establishing a practical adaptive framework for high-fidelity simulation of problems with localized sharp features in applied mechanics and computational physics.

Range reduction plays a crucial role in the accuracy and performance of evaluating trigonometric functions, and is often the primary bottleneck for large floating-point inputs. While fast algorithms such as Cody–Waite work efficiently over narrow intervals, the Payne–Hanek algorithm remains the standard technique for accurate reduction across large floating-point inputs. However, many existing implementations of Payne–Hanek suffer from high latency due to heavy branching, conversion overheads, and the use of multi-word integer arithmetic, which hinders SIMD vectorization. In this paper, we analyze and present a branch-free variation of the Payne–Hanek algorithm using only floating-point arithmetic. Our method operates directly over large double-precision inputs ($|x| \ge 2^{16}$) and is well suited to hardware with FMA instructions. We formulate the precision constraints in terms of a truncation error budget, construct a compact lookup table indexed by the input exponent, and prove that the scaled reduced argument has absolute error below $2^{-110}$ and relative error below $2^{-60}$ for every input, including worst cases. The same routine can serve both as the complete range reduction of a single-stage implementation and as the fast path of a correctly rounded one, achieving higher throughput than existing implementations and lower latency than those returning a double-double reduced argument. The algorithm is currently implemented in the LLVM libc project.

For functions that are convex of higher order, classical extremalities identify Gaussian quadrature as a lower bound and Gauss–Lobatto quadrature as an upper bound for the integral. We sharpen this bracket in every odd order: if $f$ is $(2n-1)$-convex on $[-1,1]$, then the exact integral $I[f]:=\int_{-1}^{1}f(x)\,dx$ lies not merely between the $n$-point Gauss–Legendre rule $G_n[f]$ and the $(n+1)$-point Gauss–Lobatto rule $L_{n+1}[f]$, but already between $G_n[f]$ and their midpoint. This yields the certified estimator $Q_n[f]:=\frac{3}{4}G_n[f]+\frac{1}{4}L_{n+1}[f]$, with $|I[f]-Q_n[f]|\leqslant\frac{1}{4}|L_{n+1}[f]-G_n[f]|$, and the constant $\frac{1}{4}$ is sharp. We also record a complementary phenomenon in even order: for $2n$-convex functions the Gauss–Radau endpoint rules still bracket the integral, but the midpoint of the two Radau rules cannot yield a one-sided refinement. The obstruction is explained by an odd, sign-changing Radau quadrature kernel. Finally, we distinguish numerical certification from ordinary high-order approximation: we include concrete examples from spline, moment and statistical models and report experiments for truncated powers and near-pole Stieltjes kernels, where the shape certificate remains meaningful when classical smoothness-based error constants are unavailable or severely pessimistic.

This work deals with the numerical solution of parametric linear elliptic partial differential equations using the finite element method for spatial discretization and a proper orthogonal decomposition or Lagrange subspace method based ROM for treating the parameter dependency. We study a simple one-parameter model problem and prove an error estimate for these ROMs as a function of the snapshot point set. Sub-exponential convergence with respect to the number of snapshot points is established under very mild assumptions on the point set. Exponential convergence is proven for structured point sets. In contrast to existing error analysis of Lagrange subspace methods, our estimates are not based on interpolation or perturbation arguments. Rather, we define a new approximation operator from the snapshot subspace and prove error and stability estimates for it.

Data analytics offer a great opportunity for organizations within the cold-chain industry to become more sustainable and reduce their environmental footprint. Yet, despite recent progress in data sciences, they still face challenges to ensure data quality and streamline data processes through adapted full-scale ecosystems, limiting their ability to perform sustainable analytics. Most data techniques and tools stem from the data-science community and are not always suitably adapted to domain specificities. For instance, nearly all data-quality frameworks are designed to individually consider entities, i.e. representations of real-world elements. However, achieving environmental impact reduction in industrial settings, particularly within the cold-chain sector, often requires modelling and studying entire systems comprised of multiple interdependent entities. This publication challenges the application of conventional data science and engineering concepts in data ecosystems within the cold-chain sector and proposes an alternative approach to sustainably address prevalent challenges and thus optimize analytical processes.

We present new fully discrete Cartesian grid Active Flux methods for the two-dimensional Euler equations of gas dynamics. A key component of these methods is the evolution operator used to evolve point-value degrees of freedom. Therefore, we derive several new, third-order accurate, multidimensional evolution operators for point-value updates of the linearized Euler equations using the method of bicharacteristics. We also investigate how the reconstruction and choice of approximate evolution operator affect the accuracy and stability of the resulting methods.

In this work, different from previous results, the explicit expression of numerical Onsager-Machlup action functional for Euler discretized SDEs driven by fractional Brownian motion is derived provided the drift coefficient and numerical reference path satisfy some suitable conditions. Then numerical fractional Euler-Lagrange equations for numerical Onsager-Machlup action functional are also obtained. Finally, numerical experiments to illustrate and support the theoretical findings.

We develop and analyse a symmetric weighted interior penalty polytopal discontinuous Galerkin (SWIP-PolyDG) method for multi-species reaction–diffusion systems with nonlinear reaction terms of up to cubic order. Such terms arise naturally when higher-order interactions are incorporated into population models, allowing non-additive effects among multiple species to influence local growth and conversion mechanisms. The proposed framework accommodates heterogeneous and possibly anisotropic diffusion tensors on general polygonal meshes. To control the nonlinear coupling, we consider a hierarchical block structure in the reaction operator, whereby each population block depends only on its own variables and on those associated with preceding blocks. In addition, the cubic self-interactions within each block are assumed to have a dissipative diagonal structure. Under these hypotheses and proceeding recursively over the hierarchy of population blocks, we derive local-in-time stability estimates in two spatial dimensions for the semi-discrete formulation in both the $L^2(\Omega)$-and dG-norms. We further establish an a priori error estimate in a combined $L^2(\Omega)$-dG energy norm for sufficiently regular solutions for the semi-discrete formulation. Finally, numerical experiments confirm the predicted convergence behaviour and illustrate the robustness of the method under heterogeneous diffusion.

We lay out the analysis of the “Kronecker preconditioners” adopted in [24] for the linear wave equations. Importantly, we deliver the analysis justifying the claims in~[24, Section 4.2] about these preconditioners and we unpack all of the necessary numerical background for a full understanding of the observed phenomena.

We focus on systems of nonlinear equations, where some equations are non-polynomial. Building on the remarkable developments in solving polynomial systems during the recent decades, and on a few extensions to specific examples of the transcendental case, a general framework for finding isolated roots is proposed. The name paxotopy is proposed after (1) Polynomial ApproXimation of the transcendental functions, then (2) solving the approximating polynomial system by the homOTOPY continuation method, and finally (3) the roots of the polynomial systems are used as starting points of a Newton's iteration to arrive at the solutions of the original system. These steps combine the advantages of local and global methods.

BF16 is now standard in large-scale pretraining, including in fused attention kernels such as FlashAttention, and these kernels are widely trusted. When we used FlashAttention-3 to pretrain a 450M-parameter transformer on 50B tokens, however, we ran into a problem: training was healthy for 25B tokens, then the gradient norm grew a thousandfold and the loss ended 0.2 nats above FP32 attention, without a single NaN. Recomputing the attention backward of just two layers in FP32 removes almost all of the excess gradient. Part of the cause is known: a fused multiply-add in the forward softmax, so far treated as an extreme-input NaN case and never fixed in FlashAttention-3. Repairing it stops the blow-up, but the query gradient is still wrong by more than its own size, and training still drives attention logits to thousands of times their size under accurate gradients. The remaining error comes from a broken conservation law. The softmax score gradient sums to zero along every row, which makes the query gradient blind to where the keys sit as a group; rounding it to BF16 leaves a small nonzero sum that leaks the mean key into the gradient, and the leak grows exactly as late training makes keys large and attention sharp. We introduce GProj (gauge projection), which restores the zero sum after the cast with two rank-one corrections per row. It cuts the remaining median query/key gradient errors from 219%/13% to 0.34%/0.37%, on par with FP32 attention, for 4.7% more time per training step. In matched from-scratch runs it trains to the same loss as FP32 attention, while FlashAttention-3 and key smoothing both destabilize.

Transcatheter edge-to-edge repair (TEER) treats mitral regurgitation with a clip that grasps and coapts the mitral valve (MV) leaflets. Predicting its outcome requires resolving the interaction between the blood flow and the complete MV apparatus, in which highly deformable leaflets and chordae tendineae undergo large deformation and complex contact. In this work, we introduce an immersogeometric fluid-structure interaction (FSI) formulation in which the coupling among the valvular components, the surrounding flow, and the repair is enforced weakly, without requiring conforming discretizations. A fully parameterized reconstruction pipeline is developed to convert segmented images into analysis-ready isogeometric meshes. The leaflets and chordae are modeled as Kirchhoff-Love shells and cables, and the shell formulation is modified to include an in-vivo prestrain tensor that accounts for the loading already present in the imaged geometry. We propose a subcell approach that improves the quadrature for the interfacial constraint and eliminates the over-refinement otherwise needed to prevent unphysical leakage through the leaflets. The MV repair is modeled by a penalty formulation that enforces leaflet coaptation over a prescribed region, allowing clip number and placement to be varied without remeshing or modeling the device. A lumped-parameter model at the outflow tract of the left heart imposes a flow-dependent pressure rather than a prescribed waveform, so that the split between forward and regurgitant flow matches physiological behavior. We validate the reconstruction against the segmented imaging data, demonstrate mesh convergence of the proposed formulations, and apply the framework to a prolapsed MV and its single- and double-clip repairs to quantify pre- and post-operative changes in valve dynamics and left-heart hemodynamics.

We prove polynomial-in-time Sobolev bounds on the approximation error between the twisted discrete nonlinear Schrödinger equation and its magnetic continuum limit. The analysis exploits complete integrability of the continuum flow through its Birkhoff representation to control approximation by a non-integrable Hamiltonian finite-difference lattice. We establish a focusing-defocusing dichotomy in uniform polynomial control: the bounds hold for arbitrary defocusing data and focusing data of sufficiently small mass, while modulational instability and aliasing yield an obstruction for large focusing data. The accompanying algebraic spatial rates are generically optimal in the Baire-category sense. The proof combines Birkhoff coordinates with modified energies for discrete Sobolev growth, while Fourier filtering recovers uniform spacetime estimates below the energy space despite degenerate lattice dispersion.

We present a new class of near-linear algorithms for efficiently integrating general tensor fields defined on trees with distance dependent kernels, the Structure-Adaptive Tree Field Integrators (STAD-TFIs). STAD-TFIs exploit the tree's underlying structure through decompositions built around path backbones and single vertex separators, and use two-dimensional fast Fourier transforms to compute interactions jointly. By exploiting this structural information, STAD-TFIs achieve more computationally efficient integration than their regular efficient tree field integrators (TFI) counterparts. We provide a detailed theoretical analysis of our proposed approach and complement it with an exhaustive empirical evaluation, ranging from speed tests on synthetic trees, through accelerated Sinkhorn-based relaxations of the Optimal Transport algorithms on real meshes, to Topological Attention Transformers for vision tasks. To the best of our knowledge, we provide some of the first results showing that efficient to compute and accurate relaxations of the geodesic Sinkhorn-based solutions of the Optimal Transport problem can be derived by applying fast TFI methods.

Neural operators can accelerate repeated nonlinear mechanics calculations, but their accuracy can deteriorate as operating conditions move beyond the training range. This work studies whether high-fidelity solutions acquired during use can be reused to adapt a neural operator and improve subsequent mechanics-based command selection. Two hard-magnetic soft-material systems are simulated using high-fidelity finite-element (FE) models, providing reference solutions for evaluating surrogate predictions and selected commands. A neural operator predicts deformation from known material, loading, and magnetic-field inputs, while an empirical error estimator determines which predictions may be used for command selection. Selected FE evaluations supplement these predictions, and their complete loading paths are retained for periodic updates of the neural operator and estimator. In both examples, the fixed operator loses substantial accuracy when stiffness and loading move outside the training range. Updates using 16 acquired paths recover much of the lost accuracy while preserving accuracy in the nominal regime. Under the same FE evaluation budget, the updated operators also improve command selection, although the benefit varies with the operating condition. Error estimation is less consistent, with inaccurate predictions sometimes accepted and accurate predictions rejected. These results demonstrate that reusing high-fidelity loading paths can extend the useful operating range of a neural operator. However, improved forward accuracy alone does not guarantee reliable prediction acceptance, highlighting prediction-specific error assessment as a separate requirement for trustworthy decision making.

Differential-algebraic equations (DAEs) present a ubiquitous model class in the automated modeling, simulation, control and optimization of complex finite or infinite dimensional dynamical systems, in particular in the construction of digital twins. When combined with adequate space and time discretization methods they present the work-horse of developments in many engineering domains. However, this approach comes with certain challenges for the numerical methods and in particular in the coupling of systems from different physical domains with different scales or regularity requirements. We survey the modeling with DAEs and address how to overcome some of the challenges via the framework of port-Hamiltonian differential-algebraic systems (pHDAEs) and their simulation and optimization with space-time-model adaptive methods.

Plug-and-Play (PnP) algorithms are a class of iterative methods for inverse imaging. Within an optimization algorithm, they combine a flexible fidelity term, encoding the forward operator, and a pretrained image denoiser, in order to deal with more severe corruptions such as blurring or downsampling when reconstructing an image. This work studies provably convergent PnP methods for mixed-noise forward processes by using the infimal convolution as a fidelity term, providing a statistical interpretation as a joint maximum a-posteriori estimator over both noises, and preserving the Bayesian MAP interpretation of PnP methods. Independently, we extend the PnP formulation to multiple prior terms using the Davis–Yin three-operator splitting. This extension can be combined with either standard fidelities or the proposed mixed-noise infimal-convolution fidelities. We verify that these generalized PnP methods are convergent under standard Kurdyka–Lojasiewicz conditions. Numerical experiments on Laplace-Gaussian and Poisson–Gaussian noise demonstrate stable convergence of single-prior and multiple-prior PnP methods, and divergence under fidelity mismatch. Furthermore, PnP with infimal convolution fidelities are able to scale to noise up to 38% standard deviation, with multiple priors reaching different stationary points that qualitatively preserve more textural properties.

In this work, a globally stiffly accurate Implicit-Explicit (IMEX) Runge-Kutta scheme is developed and implemented in the GBS code [Ricci et al., Plasma Phys. Control. Fusion, 2012], for two-fluid plasma turbulence simulations. The stiffest phenomena, governed by shear Alfvén waves and parallel diffusion, are treated implicitly, while the remaining non-stiff terms are advanced explicitly. This splitting enables time steps well beyond the Courant-Friedrichs-Lewy limit, without incurring the full computational cost of a globally implicit formulation. To efficiently solve the implicit subsystem at each time step, a three-dimensional physics-based preconditioner, inspired by techniques developed in the magnetohydrodynamic (MHD) context, is used. The resulting framework is verified through the method of manufactured solutions and exhibits both algorithmic and parallel scalability. Significant advantages in numerical stability and computational efficiency are demonstrated with respect to an adaptive explicit Runge-Kutta scheme.

Smoothing empirical data can fill gaps between observations, while noise in the surrounding space can move mass away from the manifold supporting the data. We derive Wasserstein error bounds for ambient Gaussian smoothing and for smoothing through an encoder–decoder representation. For a compact connected smooth manifold without boundary and a smooth positive target density, established heat-flow and entropy estimates imply that sufficiently small intrinsic heat smoothing improves the squared intrinsic Wasserstein error. Coupling this estimate with ambient noise yields a bound whose quadratic term includes ambient dimension and a curvature correction. Our main latent bound tracks reconstruction error, tangential distortion, and the decoder's response to orthogonal latent noise. Its quadratic surrogate gives bandwidth rules and explicit conditions under which the latent upper-bound surrogate is smaller than the ambient one. These surrogate comparisons do not assert an ordering of the exact Wasserstein errors. Synthetic experiments examine the predicted scaling and latent geometry on spheres and a product of spheres.

Automatic relevance determination (ARD), the default tool for variable selection in Gaussian-process (GP) regression, ranks inputs by inverse lengthscales – which measure how fast a function varies, not how much an input contributes to prediction – and offers no calibrated rule for deciding which inputs to keep. The prediction-centred alternative, the derivative sensitivity $\nu_j = \mathbb{E}[(\partial f/\partial x_j)^2]$, is available in closed form from a fitted GP, but turning it into a selection rule is harder than it looks: at a null input the estimator is a degenerate quadratic form, so Wald and Bernstein-von Mises cutoffs are anti-conservative, and the natural residual bootstrap is mis-scaled. We show that a studentized multiplier bootstrap of the GP derivative process repairs both, prove its validity through an invariance principle for quadratic forms, and obtain asymptotic family-wise and false-discovery-rate control across inputs. Over 100 replications the rule controls FDR wherever inputs are truly null, while uncalibrated derivative rankings breach the target by up to 2x and a Bernstein-von Mises cutoff by 2.2x; at matched FDR it loses no power; it holds under a Matérn kernel and input correlation up to 0.99; on real data with planted and authentic null inputs it admits 5-12x fewer spurious inputs; it costs 5-18% of the GP fit; and a block-averaged variant retains validity at cost linear in $n$.

In this work, we propose several novel numerical schemes that synergistically integrates the $L1$ method and symmetric fractional-order reduction (SFOR) method, specifically designed for fractional sub-diffusion equations with time delay. Crucially, we rigorously establish an optimal convergence order of $2 - \alpha/2$ in temporal direction, which improves upon the classical rate of $2-\alpha$ under identical smoothness conditions. To further enhance computational efficiency, we develop a weighted alternating direction implicit (ADI) scheme for the two-dimensional case and analyze its coupling with compact spatial discretization—thereby achieving high-order accuracy in space. Moreover, we introduce a new framework for analyzing the error of $L1$-type discretization on graded temporal meshes, which relaxes the mesh grading requirement and still guarantees optimal convergence rate at continuous time points. Some numerical experiments comprehensively validate both the theoretical accuracy and computational efficiency of the resulting schemes.

Deterministic neural networks and neural operators provide point predictions with no intrinsic measure of reliability. Yet, predictive uncertainty may stem from irreducible outcome variability, finite data, or limitations of the chosen model class. Monte Carlo dropout offers a computationally convenient way to construct a predictive distribution through stochastic feature masking, without training multiple independent networks or explicitly inferring a posterior over model parameters. However, its perturbation law is largely prescribed a priori and typically factorised across latent coordinates. We introduce Optimal Transport Dropout (OTD), which instead learns the predictive mapping and the law of its latent perturbations jointly. Starting from a simple independent reference distribution, OTD transports latent perturbations through a learnable flow and propagates them through the predictive neural network, thereby inducing a structured predictive law. Training uses the strictly proper Energy Score, while a kinetic-action term geometrically regularises the transport. Synthetic benchmarks show that OTD captures multimodal predictive distributions, generates meaningful dispersion when the model is misspecified, and exhibits contracting dispersion as more training data or greater model capacity are provided. For a field-valued partial differential equation surrogate, predictive dispersion strongly aligns with the spatial pattern of prediction errors. On this task, compared with Monte Carlo dropout, OTD yields more accurate predictions and better-calibrated, substantially narrower intervals. On real-world regression benchmarks, it further shows competitive accuracy and better probabilistic predictions compared to established baselines. OTD therefore offers a way to learn structured predictive uncertainty without explicit posterior inference or ensembles of independently trained predictors.

We propose an approach using matrix symmetry to accelerate hierarchical low-rank compression for wave transmission scattering. For Helmholtz transmission problems, we introduce a symmetric Müller-like boundary integral equation formulation. Although not of the second kind, this formulation creates no difficulties for direct solvers. A one-point corrected quadrature in a Nyström discretization preserves the complex symmetry of the discretized system. This formulation is integrated into a variant of a hierarchical semiseparable direct solver. Within this framework, the computational workload of the low-rank approximation in the proxy method is halved by exploiting the symmetric relationship between the left and right coefficients. Furthermore, we extend this symmetry-exploiting approach to elastic transmission problems by using a symmetric Galerkin boundary element method. Numerical experiments on a multi-node computing cluster demonstrate that the symmetric Müller-like formulation for Helmholtz scattering achieves on average a 1.41-fold speedup compared to the non-symmetric baseline while maintaining third-order convergence. Furthermore, the symmetry-aware solver for elastic scattering is more than 1.5 times faster.

Pricing European options under local volatility requires repeatedly solving a PDE whose coefficients change with recalibration, while practitioners need both prices and sensitivities across spot-time surfaces. Neural surrogates can amortize these solves, but near expiry the solution loses regularity, making curvature difficult to learn and allowing violations of no-arbitrage bounds. We propose a finance-informed Deep Operator Network (FI-DeepONet) that decomposes the pricing operator into an input-dependent Black-Scholes carrier and a learned correction. The carrier uses strike-line integrated variance, capturing the leading near-expiry curvature singularity, while a smooth monotone Financial Admissibility Layer enforces pointwise price bounds. We derive a short-maturity estimate and differentiated asymptotics for the exact correction, together with identities linking correction error to price, sensitivity, and PDE-residual errors. On in-distribution tests, FI-DeepONet reduces global relative price error by nearly an order of magnitude versus vanilla and physics-informed DeepONet baselines. On parameter-shift OOD tests within the same local-volatility family, generalization is moderate. On index-option data, the frozen model remains accurate without market-specific retraining, though it does not outperform the same-input analytic formula. We also give a quantum-compatible realization in which selected linear maps are exactly compiled or approximated by a restricted diagonal-orthogonal family before supervised adaptation.

In this paper we fully categorize all hybrid-variable (HV) discretizations of linear advection equations according to their stability property and establish the stability barrier of these methods. In the HV discretization framework, we find numerical approximations to both cell-averaged solutions and nodal solutions and evolve them in time simultaneously. We prove that the HV methods have a similar stability barrier as the classical finite difference schemes, that the semi-discretized HV scheme is stable if and only if it uses a stencil such that the number of unknowns in the upwind direction is precisely that in the downwind direction plus one or two, with only one exception when the two numbers are three and zero, respectively. It is also proved that all central HV schemes are neutrally stable, in the sense that the L2-norm of the solution remain uniformly bounded at all times. The theoretical predictions are verify by extensive numerical tests.

We study in this paper an unsupervised, two-step, model-informed deep learning framework for solving the phaseless inverse scattering problem. The objective is to reconstruct a compactly supported function that characterizes a scatterer from boundary measurements of the modulus of the total wave corresponding to multiple incident waves. In the first step, the phaseless data are transformed into an imaging function that is directly related to the underlying scatterer. Motivated by the contrast source method, we derive a coupled system of model equations using the Lippmann-Schwinger integral equation and a spectral-based equation involving the imaging function. In the second step, the unknown scatterer and contrast source function are each parameterized as independent feedforward neural networks, which are trained simultaneously at each iteration using the derived system of model equations. The Lippmann-Schwinger integral equation enforces the underlying physics as a model constraint, while the spectral-based equation incorporating the imaging function supplies geometric information on the shape and location of the unknown scatterer. The resulting framework achieves accurate and efficient reconstructions while maintaining robustness to measurement noise. Furthermore, transforming some aspects of the problem to the spectral domain enables an effective reduction in the computational cost, allowing us to recover the unknown scatterer quickly. Numerical experiments in both two- and three-dimensional settings demonstrate the effectiveness, stability, and practical potential of the proposed approach.