arXane › Computer Science › Computer Science
Wednesday, 4 November 2026: 21 new papers, 3 cross-lists and 17 replacements in Numerical Analysis (cs.NA), in arXiv's order, with full abstracts. The same papers as arXiv's own listing, set as a newspaper; every title links to arXiv.
We present a preconditioning pipeline that accelerates Tensor Train (TT) linear solvers for symmetric positive-definite problems. Extending the Factorized Sparse Approximate Inverse framework to tensor networks, we construct a preconditioner as the product of a TT operator and its transpose. This formulation is positive semidefinite for any bond dimension and truncation, and definite whenever the factor is non-singular, satisfying the mathematical requirements of the Conjugate Gradient method. We detail the construction of this preconditioner and its integration into Density Matrix Renormalization Group (DMRG)-like linear solvers. Performance evaluations show consistent iteration reductions for anisotropic diffusion and Poisson problems across various configurations, leading to net speedups of up to 25% in ill-conditioned settings of time stepping schemes that can reuse the approximate inverses. These benchmarks were performed with the open-source Julia library DifferenTT.jl, which we also introduce in this work to provide out-of-the-box access to these preconditioned pipelines.
In this paper, we introduce a generalized weighted randomized Kaczmarz framework to approximate solutions to overdetermined linear systems whose measurements are perturbed by sparse, additive corruptions. The proposed framework samples equations with probabilities proportional to a user-defined weighting function of their corresponding residual magnitudes. Our work focuses on weighting functions that will reduce the influence of potentially corrupted equations, and provides analyses for a general class of non-increasing weighting functions. We derive one-step expected squared error bounds that separate progress on uncorrupted equations from the contribution of corrupted equations, and provide sufficient conditions for geometric convergence under one specific family of weighting functions -- threshold weighting. We also investigate adaptive exponential weights scaled by the median or median absolute deviation of residual magnitudes, requiring no prior knowledge of the corruption fraction or magnitudes. Lastly, we consider a family of inverse residual weighting functions. Numerical experiments demonstrate the benefit of such a generalized framework in its flexibility in developing approaches for settings in which information about the corruptions is known or unknown.
Score-learning methods generally lack pointwise supervision, which affects the training accuracy, particularly in regions where SDE trajectories are rare. To address this challenge, we introduce a novel method called \emph{BEL anchors}, which uses Monte Carlo score estimates at user-selected space-time query points. By combining a Feynman--Kac representation of the Fokker--Planck equation with the Bismut--Elworthy--Li formula, our estimator requires neither closed-form transition densities, spatial discretization, nor derivatives of the initial density. We derive its asymptotic variance and use the anchors as sparse supervised targets for existing score-learning objectives. Experiments show that BEL anchors improve implicit score matching in poorly sampled regions, prevent score-PINN collapse, and enable accurate sample regeneration in several numerical examples.
Coarse finite-volume information does not automatically define an exact pathwise closure for the restricted evolution of a fine-grid realization. We study two one-step information boundaries that are distinct from recursive predictive memory. First, we show that parent averages together with complete initial parent-interface trace pairs are sufficient for one forward-Euler finite-volume transition with a deterministic two-state numerical flux, but they do not generally determine a completed second-order strong-stability-preserving Runge-Kutta, equivalently Heun, transition. An explicit scalar upwind collision shows that hidden interior data can affect later stage traces and hence the completed parent update, even when initial traces and parent averages agree. Second, we show that exact conservative one-step closure is not the same as entropy-accounting closure. A parent entropy balance obtained from a cellwise numerical entropy balance contains a subcell Jensen gap, a realized entropy-flux difference, and residual dissipation. For a declared Burgers-Rusanov conservative flux and a declared Rusanov-type numerical entropy flux, we construct a sign-indexed finite family in which all branches have identical parent averages, identical conservative interface registers, and identical next parent averages, while the entropy-interface register distinguishes every branch. Under the uniform branch law, the conservative channel leaves nonzero deterministic squared-error floors for both the raw entropy register and its cyclic divergence. The results give a concise finite-volume information-sufficiency test for multistage interface data and entropy-side accounting.
The Landau transformation xi = x/s(t) is the usual way to hold a moving boundary in place. It fails for fractional time derivatives because they have no chain rule: the transformed operator gains a nonlocal integral with a closed form only for self-similar fronts. We develop a sharp-interface finite-difference scheme for the conservation-consistent memory-flux fractional Stefan problem on a fixed grid, storing the memory limit h(x) = s^-1(x) as the time each node melts. The front is advanced using an integrated energy balance, with the singular source timed to the sub-step crossing instant. Against a known self-similar solution, the front constant is accurate to 0.02 to 0.11%, whereas competing formulations differ by 17 to 27% at alpha = 0.6. A remaining positivity violation follows from dt/t0 ~ s0^(1 - 2/alpha) and disappears on a front-locked time grid. We also derive two exact solutions with non-power-law fronts: a source-dominated solution and a Mittag-Leffler traveling wave, and verify the scheme against both. Neural-network solvers for moving boundaries are often trained in normalized coordinates; for alpha < 1, this silently changes the problem, while the change vanishes at alpha = 1, where checks are commonly performed. In the classical limit, our scheme is more accurate in 12 ms than a reimplemented moving-boundary network in 177 s.
Port-Hamiltonian (pH) systems provide a general framework for modeling a broad class of physical dynamical systems, including systems with dissipation. Constructing accurate and robust data-driven reduced-order models for general pH systems is challenging. In this work, we focus on an important subclass of pH systems satisfying reciprocity, a property arising naturally in a variety of models. Rather than seeking to preserve the full pH structure, we exploit the additional symmetry provided by reciprocity to develop an efficient data-driven reduced-order modeling approach. We show that every minimal realization of a reciprocal system is self-adjoint with respect to a generally indefinite inner product and can, up to a similarity transformation preserving the pH structure, be represented in a Krein space with a signature inner product. This observation motivates a Galerkin projection onto a tangential rational Krylov subspace with respect to the corresponding indefinite inner product. Although the resulting projection involves full-order state solutions and the unknown inner product, we show that its reduced mass and stiffness matrices can be constructed directly from transfer-function data through a Loewner formulation. We then propose a greedy data-driven procedure that adaptively selects both interpolation frequencies and tangential directions. At every iteration, the frequency is chosen according to the largest transfer-function approximation error, while the tangential direction is selected from a dominant singular direction of the corresponding matrix-valued residual. Numerical experiments for wave propagation in lossy media demonstrate improved approximation efficiency compared with block adaptive interpolation. We further apply the method to spectral mass gauging of a three-dimensional propellant tank to provide an accurate estimate of the tank volume through Weyl's law.
TriCG and TriMR are two short-recurrence iterative methods introduced by Montoison and Orban [SIAM J.\ Sci.\ Comput., 43 (2021), pp.\ A2502--A2525] for solving symmetric quasi-definite linear systems. Extensive numerical experiments have demonstrated the efficiency of these two methods, but a systematic convergence analysis is still lacking. In this paper, we derive convergence bounds for the TriCG error norm and the TriMR residual norm. Numerical experiments support the theoretical results.
We develop a hybrid tensor-train/quantized tensor-train (TT/QTT) upwind discontinuous Galerkin method for three-dimensional transient multigroup neutron transport. The central contribution is a fixed-matrix upwind construction: two spatial DG matrices serve all directions and energy groups, with signed and absolute direction-cosine factors preserving Kronecker structure. For homogeneous materials on uniform Cartesian grids, operator rank bounds are independent of mesh size, group count and direction count at fixed scattering order. We derive a conditional global step-defect bound and enforce group-integrated balance of the compressed system. Crank-Nicolson tests show second-order convergence on a separable manufactured solution whose finest final state occupies 13 kB instead of 35 TB, a compression factor of $2.7\times10^9$. A modified rank-adaptive basis-update and Galerkin integrator (TT-BUG) advances passive and active bursts with approximately half a billion unknowns. Selected runs are 15-20 times faster than AMEn at comparable verified scalar-flux errors. In the reported anisotropic active-burst configurations, TT-BUG takes 14.3 s on one processor versus 716 s for OpenSn on 24 processors: approximately 50 times less wall time and 1200 times less allocated processor-time cost. This compares reported configurations rather than matched physical errors. OpenSn provides an independent discretization check. The modified TT-BUG source-onset and saturation treatment is tested numerically; its order, convergence and general tolerance-to-error relation remain unproved.
This paper proposes and examines a completely discrete discontinuous Galerkin finite element scheme to approximate the solution of a class of strongly damped wave equations with uniform exponential decay property. In the first part, a semidiscrete discontinuous Galerkin approximation is employed in the spatial variable, while keeping time variable continuous and optimal error estimates are established while preserving the exponential decay property of the original system. Subsequently, a fully discrete scheme is developed and analyzed. Using suitable energy techniques, the consistency between the decay properties of the continuous and discrete systems is demonstrated. In particular, the constants appearing in various estimates remain uniformly bounded even as the damping parameters approach zero. Moreover, error estimates that preserve the uniform decay property are highlights of this paper. Further extension of our analysis includes a class of strongly damped fourth order evolution equations using SIPG as well as $C^0$-IP methods. Finally, several numerical experiments are performed to confirm the theoretical convergence results, verify the predicted decay behavior, and investigate the influence of model parameters on stability and computational accuracy.
We develop a numerical framework to efficiently simulate the cyclic voltammetry models based on discontinuous Galerkin (DG) spatial discretization combined with linearly implicit Euler time integration. The method is applied to two electrochemical mechanisms: a heterogeneous electron-transfer (ET) reaction and an electron-transfer reaction coupled with a catalytic chemical step (EC'). The governing time-dependent reaction-diffusion equations are expressed in dimensionless form and discretized in space using locally defined orthonormal Legendre polynomial basis functions. This discretization leads to a system of ordinary differential equations that is subsequently advanced in time using the linearly implicit Euler method. The resulting solver, denoted $\text{DG}_{\text{Impl}}$, is evaluated through numerical experiments and compared with MATLAB's $pdepe$ solver in terms of accuracy and computational efficiency. The numerical results demonstrate the convergence of $\text{DG}_{\text{Impl}}$ with respect to spatial and temporal discretizations and show that the method accurately reproduces the voltammetric responses of both ET and EC' mechanisms. For a prescribed relative error, $\text{DG}_{\text{Impl}}$ requires less computational time than $pdepe$ in the numerical experiments considered. These results demonstrate the effectiveness of the proposed DG formulation for the forward simulation of cyclic voltammetry models and provide a numerical foundation for subsequent adjoint-based parameter-estimation and inverse-problem formulations.
We present an algorithm for approximating a given matrix polynomial of degree $d$ by a matrix polynomial of prescribed rank and degree at most $d$. The method builds on recent advances in the theory of generic eigenstructures and factorizations of matrix polynomials with bounded rank and degree. We also give a specialized treatment of the case of matrix pencils (degree one matrix polynomials), motivated by its particular theoretical and practical importance. Lastly, numerical experiments are reported to assess the performance of the proposed approach and to compare it with existing methods.
We present a theoretical and numerical analysis of novel IMEX Active Flux schemes for a hyperbolic approximation of the Korteweg-de Vries (KdV) equation. For this approximation, it is important that the numerical schemes recover the limit equations (asymptotic preserving) and retain their order of accuracy in the limit (asymptotically accurate). Using formal asymptotic analysis and numerical experiments, we show that an Active Flux spatial discretisation (closely related to the standard semi-discrete Active Flux scheme), combined with an implicit-explicit Runge-Kutta scheme, is asymptotic preserving but not asymptotically accurate. We therefore propose alternative spatial discretisations that differ in the update equations of the point values. First, we derive a scheme by rescaling and eliminating terms in the upwind matrices. Second, we develop a second-order scheme by using central discretisation for the terms treated implicitly in time and present two modifications of it to third-order. Finally, we formulate a scheme with predefined upwind directions that are not chosen based on characteristics. We show that the alternative schemes are not only asymptotic preserving but also asymptotically accurate, and compare their accuracy, stability, and energy conservation numerically.
Separable nonlinear least squares problems arise in inverse problems such as semi-blind image deblurring, where the image enters the forward model linearly, and the blur parameters enter nonlinearly. Even when the image is regularized, joint estimation of the image and blur parameters can favor a degenerate no-blur solution, in which the blur operator tends to the identity and the reconstruction merely denoises the data without removing the blur. We address this by regularizing the blur parameters. Motivated by the determinant term arising from Gaussian marginalization, we construct a family of blur-parameter priors whose weight and auxiliary operators can be chosen independently of the image prior. The resulting log-determinant penalty is incorporated into regularized variable projection, which eliminates the image variables by minimization. Under periodic boundary conditions and assuming the matrices involved admit a joint spectral decomposition, we show that a sufficiently large determinant weight rules out the no-blur solution as a local minimizer. We also characterize the regularizer's behavior as the auxiliary regularization scale tends to zero or infinity. We solve the reduced problem with a quasi-Newton method that evaluates the objective and its derivatives in the Fourier domain. In experiments with Gaussian, defocus, and motion blur, the proposed prior improves image and blur recovery over Gaussian marginalization across several settings, with determinant weights computed from the true image. A single-sample experiment illustrates how to automatically choose parameters from the data using a conditional generalized cross-validation criterion.
We generalize the compact residual-screened local reduced-order model (CL-ROM) from Brinkman-penalized flow to shape optimization with the diffuse domain method (DDM). As with Brinkman penalization, DDM extends evolving-domain problems to a fixed rectangular domain, allowing one grid for discretization and directly compatible snapshots from different designs. Unlike Brinkman penalization, DDM uses a phase field to weight operators and sources and treats Dirichlet, Neumann and Robin conditions through diffuse boundary terms. The CL-ROM mechanism carries over: recent full-order and corrected solutions determine local reduced spaces, while residual screening, correction with retained factors and objective verification decide when a full-order solve is needed. We derive reduced operators, residuals and sensitivities for each diffuse model, using Galerkin projection for the symmetric problems and least-squares Petrov--Galerkin projection for nonlinear flow. For the symmetric problems we prove a local approximation estimate that admits inexact, corrected snapshots. In numerical tests of shape identification, drag minimization in Navier--Stokes channel flow and compliance minimization of a cantilever, CL-ROM uses reduced solutions for nearly all design updates, reproduces the full-order objectives with relative differences below $5\times10^{-6}$, and is $3.85$ to $7.65$ times faster than a full-order baseline that reuses factorizations as preconditioners.
Accurate streamflow estimation requires a consistent representation of model and observation uncertainty. Such consistency is difficult to obtain when rainfall information is unavailable, or not temporally aligned with the discharge observations. We develop an adaptive Kalman filter for Generalized Shot Noise streamflow models, with recursive estimation of the measurement and process-noise covariances. The process-noise update includes a correction for the covariance propagated through the state dynamics, without which the recursion has no finite stationary point. With the correction, the expected increment vanishes when the mean normalized innovation squared reaches one. The filter is assessed using 37 years of daily inflow observations from the Camastra Basin in southern Italy, under two rainfall configurations. The first supplies a synthetic effective-rainfall realization as a deterministic input, whereas the second represents the unresolved forcing through the adaptive process noise. The second configuration yields the more accurate estimates, since the synthetic rainfall events, although statistically calibrated, do not coincide in time with those that generated the observed discharge. Over the full record, it attains a Nash--Sutcliffe efficiency of 0.95 and a median relative error of about 1%. Regime-conditional consistency diagnostics are reported and discussed for the second configuration.
Summation-by-parts (SBP) operators in combination with simultaneous-approximation-terms (SATs) for boundary conditions provides a recipe that can be applied to well posed linear initial boundary value problems.The SBP-SAT frame work guarantees semi-discrete energy stability and can be applied to most numerical discretization methods. High order accuracy comes without additional difficulty. The extension of SBP-SAT to nonlinear initial boundary value problems has so far been somewhat unclear. However, recent developments have shown how to construct stable nonlinear SAT terms. In this note we develop nonlinear SBP operators. These combined developments provide an SBP-SAT based recipe for stable approximations of nonlinear skew-symmetric hyperbolic initial boundary value problems.
We develop a unified isoparametric moving mesh partial differential equation (MMPDE) method for $P_1$ simplicial and $Q_1$ tensor-product elements, including elements on embedded surfaces. Mesh movement follows a prescribed symmetric positive-definite metric through equidistribution and alignment. For $Q_1$ elements, we prove that bounded, exactly integrated adaptation energy prevents local Jacobian rank loss on a fixed finite mesh at sharp critical exponents. The result also gives local full-rank preservation for embedded bilinear surface elements. A complete physical-coordinate gradient accounts for metric variation at moving quadrature points. For the discrete MMPDE, we use whole-element validity tests and develop sufficient Bernstein certificates for trilinear hexahedra. An explicit positive-margin bound on each fixed stage is provided by an optional certificate barrier. Numerical examples illustrate metric-driven movement on one-dimensional, planar, volume, and graph-surface meshes. Limitations of pointwise hexahedral Jacobian sampling are also presented.
Model reduction via the Mori-Zwanzig (MZ) formalism is well understood for coarse-grained representations constructed from standard basis functions such as Fourier functions. We show that the same formalism can be used to build reduced-order models in bases discovered by physics-informed machine learning. The fluctuation-dissipation theorem (FDT) is used to close the resulting non-Markovian evolution equations for the expansion coefficients of the learned basis. We further introduce a complementary, generative closure model. Relative to the Markovian dynamics alone, both the MZ and generative closures substantially improve predictions of the filtered states.
In this paper, we consider three-dimensional exterior wave propagation problems with non-homogeneous Neumann boundary conditions. We reformulate the problem as a space-time boundary integral equation involving a hypersingular operator. To solve the resulting formulation, we combine Lubich's convolution quadrature with a shifted collocation boundary element discretization and investigate numerically the efficiency, stability, and convergence of the resulting numerical scheme. Numerical tests show that the proposed approach achieves an accuracy comparable to that of Lubich's convolution quadrature combined with a Galerkin boundary element method, and provide numerical evidence of the unconditional stability of the collocation scheme.
This work presents an optimization-based remap approach for arbitrary Lagrangian-Eulerian (ALE) hydrodynamics that preserves sharp features and physical bounds while conserving global quantities such as volume, mass, momentum, and total energy. Achieving these properties is the main challenge in remap methods. The presented approach is well-suited for high-order finite element (FE) discretizations and ALE methods that allow for large, purely Lagrangian displacements between successive remeshing and remapping steps. Recent advances in parallel, GPU-capable interpolation routines enable the initial transfer of all fields via direct physical-space interpolation, avoiding computationally costly and diffusive pseudo-transient PDE-based advection. The interpolation produces bounded fields with minimal numerical diffusion. However, it may violate the desired conservation of global properties. Thus, a subsequent optimization step restores conservation while preserving physical bounds. The presented method provides minimally diffusive remapping of finite element fields of arbitrary order. Moreover, it directly remaps quantities defined at integration points, eliminating the need for projection onto finite element spaces. The approach also mitigates typical challenges associated with synchronizing primal variables and conserved quantities in advection-based remap methods. The proposed optimization-based remap method is evaluated using standard two- and three-dimensional remap benchmarks and full ALE hydrodynamics simulations.
Single precision (FP32) halves the memory and data movement of double precision (FP64), and on accelerators is often several times faster. Its rounding errors are harmless except in the stored state. A modified equation for Runge-Kutta methods in Shu-Osher form and linear multistep methods separates the rounding of the stored state, a forcing proportional to $1/\Delta t$, from truncation and from right-hand-side rounding, which sets a $\Delta t$-independent floor. We bound this forcing for strong-stability-preserving schemes and identify two mechanisms that keep its sign fixed, so it accumulates as a bias: past a turnaround step, refining the step increases the error, in FP64 as in FP32. For a model oscillator the analysis predicts the floor and the FP32-to-FP64 ratios of turnaround step and minimum error within 30%. Our remedy needs no FP64 arithmetic: each state variable is stored as the unevaluated sum of two FP32 numbers, updated with error-free transformations, while every right-hand-side evaluation stays in FP32. Unlike compensated summation, it needs no rewriting of the scheme. In model problems from geophysical flow, mechanical dynamics and celestial mechanics, with five Runge-Kutta and multistep schemes, it restores design-order convergence to the floor. In large-eddy simulations of a high-lift airplane wing on three refined grids, up to 544 million cells, with time steps $\Delta t$, $\Delta t/2$ and $\Delta t/4$, FP32 overpredicts the lift by 2.1, 4.5 and 7.7%, while the two-float lift agrees with FP64 within 0.3% on every grid, and on graphics processing units the production solver runs 1.6 to 2.7 times faster than FP64.
We study the multi-modal Laplace approximation, given by a weighted Gaussian mixture centred at the local minima of the negative log-density, for sampling from probability measures that concentrate around finitely many modes as an inverse-temperature parameter $n$ grows. Such distributions occur in Bayesian inference in the small noise or large data regime. We derive explicit, non-asymptotic bounds on the Hellinger distance between the multi-modal Laplace approximation and the concentrating posterior. The bound decays like $n^{-1/2}$ with $n$ and the constants depend only polynomially on the dimension $d$. Moreover, we analyse the Laplace mixture as a proposal distribution for self-normalised importance sampling and independent Metropolis--Hastings (IMH) targeting concentration multimodal distributions. We derive asymptotic results, as $n\to\infty$, on the effective sample size in case of importance sampling as well as on the average acceptance probability, lag-one autocorrelation, and spectral gap in case of IMH. Numerical experiments illustrate and verify our theoretical findings.
We study the obstacle problem for the evolution equation $\partial_t^\alpha u+(-\Delta)^s u=0$ in the whole space $\mathbb{R}^d$, where $\partial_t^\alpha$ with $\alpha\in(0,1)$ is the Caputo derivative in time and $(-\Delta)^s$ with $s\in(0,1]$ is the fractional Laplacian in space ($s<1$) or the classical Laplacian ($s=1$). We prove that the problem has at most one bounded viscosity solution, and we obtain existence as the limit of an explicit monotone finite-difference scheme, which combines the L1 discretization of the Caputo derivative with powers of the discrete Laplacian. The scheme is also used to establish regularity estimates in space and time. We then study the large-time behaviour: solutions converge, along subsequences, to stationary solutions, and we identify the limit when the initial datum coincides with the obstacle (the minimal stationary solution) or has a limit $c_*$ at infinity. In the latter case the dimension plays a role, as a consequence of the recurrence of the underlying process: the limit is the unique stationary solution with limit $c_*$ at infinity if $d>2s$, and a constant (the maximum between $c_*$ and the supremum of the obstacle) if $d\le2s$. Numerical experiments in one dimension, based on an explicit stationary solution, illustrate the convergence of the scheme and the relaxation towards the stationary profile.
We attempt to alleviate failure modes of classical physics-informed neural networks (PINNs) in scenarios of topological pathologies corresponding to punctures, holes, poles, and voids by introducing de Rham cohomological augmentations, or Hodge-de Rham PINNs. Vanilla PINNs often encounter pitfalls in hole-based settings due to discontinuous interfaces, nonzero closed path integration, path dependence, correspondence to exact forms via the exterior derivative, etc., therefore we absolve traditional failures via parameterizing a 1-form by the Hodge-Morrey-Friedrichs (HMF) decomposition for manifolds with boundary, drawing a connection to de Rham cohomology such as through the HMF isomorphism theorem. 1-forms can act as representations of scalar-valued PDE solutions through integral operator mappings. We propose performing a PINN-based loss on the parameterized 1-form through an associated PDE, which subsequently solves the potential-based PINN problem and moreover remedies the failures of vanilla PINNs via the Sobolev-regular operator under a Fr\'echet differentiation interconnection. The operator relationship often has a nontrivial closed-form and often requires integration or a case-by-case analysis, which we address thoroughly in our experiments. When an integrability condition is not met on our losses, we can equip exhaustion functions to train the networks. We demonstrate our methods on the Aharonov-Bohm effect, which, to our understanding, is largely nonexistent in PINN literature. We develop theoretical guarantees of our methods and demonstrate our pathological safeguards on a broad selection of topologically-adverse problems, including scenarios of inexact forms, and those of nontrivial harmonics and coexact terms.