Peer Revenue Fully funded

Papers

Fresh from arXiv, typeset properly. Highlighted words trade on the market.

Mathematics Numerical Analysis

math.NA Watch Numerical Analysis (sign in to use)

Thursday 13 August

Numerical Analysis

Finite element approximation of the linearized Cahn-Hilliard-Cook equation

The linearized Cahn-Hilliard-Cook equation is discretized in the spatial variables by a standard finite element method. Strong convergence estimates are proved under suitable assumptions on the covariance operator of the Wiener process, which is driving the equation. The backward Euler time stepping is also studied. The analysis is set in a framework based on analytic semigroups. The main effort is spent on proving detailed error bounds for the corresponding deterministic Cahn-Hilliard equation. The results should be interpreted as results on approximation of the stochastic convolution, which is a part of the mild solution of the nonlinear Cahn-Hilliard-Cook equation.
Read the full abstract

Numerical Analysis

BTF-PINN: Enforcing Dirichlet Boundary Conditions Without Boundary Training

The homogeneous Dirichlet boundary value problem captures the core difficulty of solving Dirichlet problems with non-interpolatory methods. We propose BTF-PINN (Boundary-Training-Free Physics-Informed Neural Network), an interior-only strategy for solving homogeneous Dirichlet boundary value problems that requires no boundary training, boundary penalties, or boundary-conforming parametrizations. The key idea is to embed the essential boundary condition into a newly designed boundary-free loss function. We prove the equivalence between the proposed interior-only variational formulation and the original boundary value problem, establish a sharp threshold condition for the residual weight, and develop a convergence analysis based on the coercivity of the functional. Numerical experiments on high-dimensional problems, irregular geometries, and anisotropic elliptic equations demonstrate the effectiveness of BTF-PINN. Comparisons with standard boundary-penalty PINNs further show that BTF-PINN achieves superior boundary trace accuracy.
Read the full abstract

Numerical Analysis · Cross-listed from math.NT

Algebraic constructions of point sequences with quasi-uniform two-dimensional projections

Motivated by sequential space-filling designs for computer experiments, we study algebraic constructions of extensible point sets in the $d$-dimensional unit cube whose two-dimensional coordinate projections are all quasi-uniform. Our two constructions share a common parametrization in terms of finite configurations of distinct rational directions on the projective line $\mathbb{P}^1(\mathbb{Q})$. First, using a cubic number field, we construct explicit Kronecker sequences for which the mesh ratios of all two-dimensional coordinate projections remain uniformly bounded over every initial segment of length $N\ge 2$. Second, using a real quadratic field, a split prime, and a compatible $p$-adic embedding, we construct nested rank-1 lattice designs with the same uniform projection property at every nesting level. The proofs combine algebraic norm estimates with transference principles between simultaneous and dual Diophantine approximation, yielding lower bounds for the separation radii and upper bounds for the covering radii, both of optimal order in the number of points, uniformly over all coordinate pairs. We also investigate how the choice of rational projective coefficients affects the resulting mesh ratios. This leads to a minimax problem for finite configurations on $\mathbb{P}^1(\mathbb{Q})$, in which one seeks to minimize the maximum mesh ratio over all two-dimensional coordinate projections. These constructions provide extensible point sets with uniformly controlled bivariate geometry.
Read the full abstract

Wednesday 12 August

Numerical Analysis

Identifying changing partial differential equations using Sampled Local WeakIdent

We propose Sampled Local WeakIdent (SLW-Ident), a framework for identifying changing governing equations from a single set of given data. Different from a typical approach of using finite element based approximation to represent varying coefficients, this paper explores a local approach in identification of differential equations. First, we present the power of Local WeakIdent which gives good local identification and is also computationally efficient with a small patch size, yet it can be sensitive to local perturbations. We propose SLW-Ident which stabilizes the identification process and also incorporates global information: we first sample patches in the whole given domain, identify equations for each sampled patch, then use residual error of these equations to find the transitions between different equations. We refer to a region where the support of the identified equation does not change to be a region of one equation. Within each region of one equation, we pick the most frequently identified equation as the identified equation, and find constant as well as varying coefficient PDEs within each region of one equation. This is justified by an uncertainty quantification theory that gives the relation between statistical error of dominant support selection and the number of patches. We provide various numerical experiments showing that SLW-Ident accurately recovers the regions of one equation and the governing equations for changing PDEs with varying coefficients even with noisy given data.
Read the full abstract

Numerical Analysis

Halving the size of skew-symmetric eigenvalue problems via the polar decomposition

This paper introduces a novel algorithm for computing eigenvalues and eigenvectors of a dense real skew-symmetric matrix $A$. Its main ingredient is the computation of a polar factor of $A$ that is both skew-symmetric and orthogonal. This polar factor is then used to transform the original problem into a Hermitian eigenvalue problem of half the size, which can be solved accurately and efficiently with standard software such as LAPACK. Numerical experiments demonstrate the stability of the method and show that its running time is competitive with existing approaches for skew-symmetric eigenvalue problems. Finally, we show that the same principle can be used to reduce an orthogonal eigenvalue problem to a unitary eigenvalue problem of half the size.
Read the full abstract

Numerical Analysis

Numerical methods for the simulation of quantum walks and quantum annealing

It is known that Chebyshev based polynomial approximation gives a near-optimal rate of convergence for calculating a function of a Hermitian matrix. However, previous work has not discussed the option of true minimax approximation, nor the specifics of writing a high performance implementation with a rigorous analysis of errors. This work provides such an analysis and an open-source implementation of three approximation methods in C++.
Read the full abstract

Numerical Analysis · Cross-listed from math-ph

Model-Consistent Structured-Hankel Matrix Completion for 3D Sparse Multi-frequency Electromagnetic Source Reconstruction

Reconstruction of 3D electromagnetic currents from sparse, multi-frequency, far-field radiation data is severely ill-posed because sparse sampling masks specific current Fourier modes and the missing spectrum enlarges the null space of the forward operator. Although structured-Hankel completion can recover missing spectrum by exploiting the finite rate of innovations structure of a compactly supported geometrically sparse source, it treats current components independently without guaranteeing compatibility with the model. This often renders non-physical solutions that violate model consistency. To resolve this, we propose a training-free Maxwell-model-consistent 3D structured-Hankel framework (MC-Hankel) that alternates joint low-rank Hankel completion with an exact, closed-form projection onto the true Maxwell synthesis subspace while enforcing conjugate symmetry and noise-aware data fidelity. Across tests on two source models, three sparse sampling rates (30%, 40%, and 50% of a Nyquist-sampled grid), and clean versus noisy (10 dB SNR additive white Gaussian noise) conditions, MC-Hankel consistently outperforms sub-sampled Fourier inversion, scale-normalized $\ell_1$-compressed sensing baseline, and the standard joint 3D annihilating filter-based low-rank Hankel matrix completion approach (ALOHA). Across all considered source-noise-sampling configurations, MC-Hankel improves mean peak SNR (PSNR) by 0.65-4.22 dB and tri-planer structural similarity (3D SSIM) by $0.0266$-$0.0703$, reduces relative full-volume $\ell_2$ reconstruction error by $7.4\%$-$38.9\%$, and brings non-physical model residuals down to machine precision over standard joint 3D ALOHA with merely 17% computational overhead. Exact signed-rank tests are underpowered for five paired trials, so statistical conclusions are reported together with paired bootstrap intervals, effect sizes, and small-sample limitations.
Read the full abstract

Numerical Analysis

A Positivity-Preserving Expectation Scheme for Hamilton–Jacobi–Bellman Equations with Oblique Robin Boundary Conditions

For anisotropic diffusion with mixed derivatives, standard compact coordinate-aligned stencils and some local finite-volume or finite-element constructions can lose nonnegative coefficients unless suitable coefficient or mesh conditions are imposed. Here nonnegative coefficients are generated directly from conditional expectation rather than from an algebraic stencil decomposition. We construct a positivity-preserving scheme for possibly degenerate Hamilton–Jacobi–Bellman equations with controlled oblique Robin boundary conditions. At interior nodes the construction stems from a one-step reflected Feynman–Kac identity: an $m$-dimensional Rademacher vector generates $P=2^m$ equally probable weak-Euler branches, an exterior branch is mirrored about its oblique projection, and twice its overshoot serves as the discrete boundary local time $D$. The Robin coefficient enters only through the attenuation factors $e^{-\kappa D}$ and $De^{-\kappa D/2}g$, leaving the equal branch probabilities unchanged. Boundary nodes use a separate same-level closure over a spatial offset $\ell_h\asymp h$, which may be implicit; uniform obliqueness and $\mathbb P_1$ interpolation give a mesh-independent positive weight on interior nodes, so fixed linear Robin data yield a sparse nonsingular $M$-matrix system and the general control set yields a monotone contraction. Positivity is preserved for nonnegative data without a diagonal-dominance condition and without any CFL-type relation between $\Delta t$ and $h$.
Read the full abstract

Numerical Analysis

Positive quadrature and mobile sampling of multivariate trigonometric polynomials

The geometric properties of quadrature points in one and multiple dimensions are a classical topic in numerical analysis. Recently, generalized quadrature methods, where discrete points and weights are replaced by integration along curves, have attracted growing interest. This paper focuses on positive generalized quadratures that are exact for multivariate trigonometric polynomials. We derive upper bounds on the covering radius of such quadratures using sign-localized test functions, a technique that also extends naturally to algebraic polynomials on intervals, spherical polynomials, and hyperbolic cross trigonometric polynomials. For lower bounds on the length of quadrature curves, we compare two recent approaches in the context of multivariate trigonometric polynomials. The second part of the paper studies a well-known, geometrically simple curve that underlies rank-1 Korobov lattice rules for periodic functions. We prove that these curves are quasi-optimal with respect to multiple optimality criteria and demonstrate its adaptability to the 2-sphere.
Read the full abstract

Numerical Analysis

Fast eight-order Pade schemes based on Chebyshev polynomials for direct Zakharov-Shabat problem

In this work, we construct fast eighth-order Pade schemes for the direct Zakharov-Shabat scattering problem. The schemes are based on an eighth-order exponential integrator obtained from the Magnus expansion. A direct extension of the conventional fast Pade representation to the eighth-order case leads to insufficiently accurate fast variants in the considered tests. To overcome this difficulty, we reformulate the spectral dependence on a finite real interval using the Joukowski mapping and represent the local numerators and denominators in the Chebyshev polynomial basis. This makes it possible to construct the global transition matrix by a fast product tree while retaining a compact polynomial representation. Numerical experiments for chirped hyperbolic secant potentials with both signs of dispersion show that the proposed Chebyshev-based fast schemes substantially improve the accuracy of the corresponding direct fast variants and can be used for efficient computation of the continuous nonlinear spectrum.
Read the full abstract

All authors Sergey Medvedev, Igor Chekhovskoy, Irina Vaseva, Mikhail Fedoruk

Numerical Analysis · Cross-listed from math-ph

A uniform elliptic reduction, an order-matching criterion, and precision benchmarks for the strong-coupling Birman-Schwinger analysis of the lattice three-boson trimer

We present a rigorous strong-coupling Birman-Schwinger analysis of the three-boson lattice Schroedinger operator on Z^2 at the exceptional quasimomentum K = pi. First, we provide an exact closed-form benchmark for the fiber Fredholm determinant, valid for every quasimomentum K, obtained via a uniform elliptic reduction. Second, we formulate and prove a general order-matching criterion determining whether a leading-order Fredholm determinant asymptotic, with relative error O(1/mu), suffices to fix the constant-order additive energy correction, or if the next-order refinement is required. Applying the criterion to the formal branch z = -2mu + d, we identify an algebraic crossing -2mu + 6 + 8/mu + O(mu^{-2}), but demonstrate that this crossing does not correspond to a true eigenvalue of the full Hamiltonian. The actual ground state obeys the rigorous variational bounds -3mu <= z_1^{pi,s}(mu) <= -3mu + 6, so that z_1^{pi,s}(mu) = -3mu + O(1), the same leading branch as at K = 0; direct finite-volume diagonalisation confirms the refined asymptotic -3mu + 6 + O(mu^{-1}) and spectral gap 2mu - 2 + O(mu^{-1}). We independently confirm that the known K = 0 constant C approximately 3.96458 requires no analogous refinement. Finally, we compare the asymptotic precision levels achieved across recent lattice few-body models and draw a structural parallel with the parity-based classification of topological band insulators at time-reversal-invariant momenta.
Read the full abstract

Numerical Analysis

Schur-Riesz Variational Enrichment: A Generalized Refinement Framework for Finite Elements

Width theory identifies economical spaces for compact PDE solution families, but does not provide a stable adaptive selection rule. We introduce Schur-Riesz refinement, which compares ordinary h/p refinement with operator-informed functions on a common variational scale. Projection removes content already represented by the incumbent, Riesz bounds test coefficient stability, and exact Schur gain ranks the surviving directions. New non-regression, bulk-contraction, and mixed near-oracle results provide finite-run error, dimension, and work certificates for coercive Galerkin and noncoercive minimum-residual formulations. We applied Schur-Riesz refinement across coercive and noncoercive PDEs. At matched dimension, automatic modes reduce held-out error by $50.2\%$ for 2D heterogeneous Helmholtz and $30.9\%$ for 2D Darcy. At matched Darcy error, the selected space uses $38.2\%$ fewer coordinates, reducing deployment memory by $38.9\%$ and online time by $30.1\%$. The additional offline construction is reusable and can be amortized across multiple right-hand sides for a fixed operator. On locked 3D heterogeneous Helmholtz tests, the method reduces mean error by $25.1$-$47.3\%$ against matched spectral-polynomial spaces and meets target errors with 24-96 coordinates, versus 1,331-6,859 for native MFEM hp-refinement. In an adaptive audit, 48 coordinates attain mean error $.1092$, compared with $.0967$ using 878 polynomial coordinates, while making online solves $6.8\times$ faster. Thus, the experiments demonstrate that Schur-Riesz refinement can construct smaller stable spaces with lower deployment memory and repeated-solve cost, while preserving a certified incumbent when richer functions do not help.
Read the full abstract

Numerical Analysis

An Interface Green's Function Framework for Complete Discrete $W^{1,\infty}$ Analysis of Discontinuous Galerkin Methods

Pointwise error analysis of discontinuous Galerkin (DG) methods for the Poisson equation has received considerable attention during the past two decades. However, on convex polyhedral domains, existing analyses can only establish optimal error estimates in the broken $W^{1,\infty}$ seminorm for several DG methods. Since the broken $W^{1,\infty}$ seminorm does not control discontinuities across mesh interfaces, a complete discrete $W^{1,\infty}$ theory for discontinuous approximations on convex polyhedral domains has remained unavailable. In this paper, we develop an interface Green's function framework for the complete discrete $W^{1,\infty}$ analysis of DG methods on convex polyhedral domains. The proposed framework introduces new interface Green's functions that represent jumps of discontinuous approximations across mesh interfaces. Its central analytical ingredient is a new local energy estimate for these Green's functions, obtained by exploiting a cancellation between neighboring discrete delta functions. This estimate differs fundamentally from existing Green's function estimates and enables us to derive maximum-norm estimate for interface jumps without introducing additional logarithmic factors.
Read the full abstract

Numerical Analysis

Contour Hankel dynamics and indicator fields for the Riemann $\Xi$-function

We develop a moving-contour Hankel framework for encoding local zero configurations of the Riemann $\Xi$-function. Weighted contour integrals of the logarithmic derivative $\Xi'/\Xi$, expressed in a holomorphic coordinate associated with the contour, are identified with the power moments of a finite atomic measure supported at the coordinate images of the enclosed zeros. This representation yields exact zero-free and rank criteria and, for conjugation-compatible contour-coordinate pairs, an inertia formula: once the matrix order is at least the number of distinct coordinate nodes, the negative index equals the number of distinct nonreal conjugate pairs. Consequently, the Riemann hypothesis admits a local finite-dimensional Hankel-positivity formulation, although establishing this positivity independently of the zero set remains unresolved. As the contour moves, the Hankel matrix evolves by a continuous congruence flow between zero crossings and undergoes finite-rank jumps at crossing events. An isolated zero produces a signed rank-one jump, whereas a nonreal conjugate pair produces a rank-two indefinite event in a real-axis circular scan that meets the pair. Removing the continuous coordinate drift yields a piecewise-constant matrix process from which crossing coordinates, multiplicities, and zero locations can be recovered. Numerical experiments validate the contour quadrature, indicator fields, continuous flow, crossing signatures, and recovery procedure.
Read the full abstract

Numerical Analysis · Cross-listed from cs.CV

New Orthogonal Multiwavelet Filters Derived by Matrix Spectral Factorization

The paper considers the construction of two new orthogonal multiwavelets with supercompact support by using the Fast Bauer's method for matrix spectral factorization on the matrix product filter of the orthogonal CL multiwavelet filter. The new multiwavelets possess orthogonality, symmetry/antisymmetry, and one of them provides better coding and smoothness than other supercompact multiwavelets. The performance of the new multiwavelet filters in subband-based edge detection, grayscale and color image compression and 1D and 2D signal denoising is compared with the GHM, SA4, CL, Integer Haar and Alpert multifilters. The comparative analysis shows that new multiwavelets can provides better human visual measures, SSIM and MS-SSIM in image compression and denoising applications.
Read the full abstract

Tuesday 11 August

Numerical Analysis

Runge–Kutta-Aligned Oscillation Elimination: Restoring Superconvergence for Fully Discrete Shock-Capturing DG Schemes

Nonlinear stabilization is indispensable for discontinuous Galerkin (DG) discretizations of hyperbolic conservation laws, yet it typically disrupts the delicate error structure required for superconvergence analysis. Consequently, existing theory has largely been restricted to linear or semi-discrete schemes lacking oscillation control. This paper bridges this gap by proposing a Runge–Kutta (RK) aligned oscillation-eliminating (OE) DG framework that restores the superconvergence properties. By synchronizing the pseudo-time step in the OE procedure with the cumulative RK stage coefficients, we unlock a cancellation mechanism for low-order interface errors that is inaccessible to standard OEDG formulations. We rigorously prove that this fully discrete scheme achieves $(k+2)$-th order superconvergence to a tailored projection of the exact solution for linear conservation laws in both one and two dimensions, while maintaining the non-oscillatory shock-capturing capabilities of the original method. Moreover, we establish a general guiding principle for designing a class of OE-type DG schemes that exhibit such superconvergence. Key theoretical innovations include the construction of stage-aligned correction functions to compensate for nonlinear OE sources and the discovery of a two-dimensional projection operator that preserves outflow-edge averages, a property essential to close the discrete shift estimates. Numerical experiments confirm the predicted superconvergence rates and demonstrate that RK alignment preserves the parameter-free robustness of OEDG for problems with strong discontinuities.
Read the full abstract

Numerical Analysis

Stability of Finite-Batch Particle Mean-Field Variational Inference Beyond Strong Convexity

We study the implementable finite-batch particle algorithm for mean-field variational inference as a fully discrete stochastic approximation of the projected Wasserstein dynamics. The target potential is globally smooth but need not be strongly convex. The departure from contractivity is quantified by the curvature defect \[ \mathfrak d_\alpha(x,y) = \bigl[\alpha\|x-y\|^2- \langle\nabla V(x)-\nabla V(y),x-y\rangle\bigr]_+, \] which is the additive loss in the one-step Euler contraction estimate. We prove a non-asymptotic Wasserstein stability bound that separates initialization, product-empirical approximation, finite-batch drift error, time discretization, and the defects accumulated along the coupled trajectories. Under the uniform bound $\mathfrak d_\alpha\leq\beta$, the particle iterates remain within $O(\sqrt{\beta/\alpha})$ of any MFVI minimizer, up to explicit errors in the particle number, batch size, and step size. The proof uses a stationary comparison array whose population law is an MFVI minimizer but whose particle-level law is a random product empirical measure, and it controls the resulting projected-drift discrepancy explicitly. We also give coordinatewise defect estimates and structural conditions for dimension-independent projected-drift sensitivity, construct an arbitrary-dimensional smooth nonconvex benchmark with a closed-form MFVI minimizer, and explain why polynomially growing drifts require a modification of the untamed explicit scheme.
Read the full abstract

Numerical Analysis

Linear-cost Polyharmonic Spline Interpolation of Arbitrary Degree

We introduce a simple and performant approach for rapidly and accurately performing polyharmonic spline (PHS) interpolation using a combination of the fast multipole method (FMM) from computational electrostatics and the Vecchia approximation from the Gaussian process and sparse approximate inverse literatures. Using basic properties about Hadamard products and low-rank matrices, we demonstrate that an FMM with two kernels, the logarithmic and distance kernels, results in fast PHS interpolation for all orders. Furthermore, we demonstrate the exceptional performance of sparse inverse approximation methods with the Matérn covariance model for preconditioning. Combined with careful management of disallowed subspaces, we describe a procedure for obtaining prediction weights using preconditioned conjugate gradient that converges in less than $15$ iterations, even for problem sizes with over one million points. As a result, thin-plate spline interpolation—a particularly popular method that does not require parameter tuning—that matches the fully dense $\mathcal{O}(n^3)$ computation in accuracy can be done at the cost of approximately $50-60$ FMMs. A high-performance software library for odd-order PHS interpolation in two dimensions is made available as a companion to this work.
Read the full abstract

Numerical Analysis

Exact characterisation of maximum-angle conditions for spherical finite element meshes

Maximum-angle conditions are standard finite-element mesh hypotheses that permit anisotropic triangles excluded by minimum-angle or shape-regularity assumptions. For exact spherical triangles, however, the angles of the chordal affine core and the intrinsic spherical angles need not coincide, while radial geometry introduces curvature-scale distortion. We give an algebraic characterisation of the intrinsic spherical maximum-angle condition through a dimensionless quantity computed from the three vertex vectors. A uniform positive lower bound on this quantity is equivalent to a uniform spherical maximum-angle bound and requires neither spherical-angle nor spherical-area evaluation. We identify the support-plane geometry linking the chordal circumradius to radial distortion and derive sharp comparisons between the intrinsic spherical and chordal semi-regularity parameters. In particular, a locality-independent comparison holds with sharp constant $2/\sqrt3$ and a characterised equality case. An area-based parameter is shown to be smaller than the spherical semi-regularity parameter, with sharp constant one in the local flat limit. For spherical finite-element meshes, the vertex criterion implies uniform chordal semi-regularity with an explicit constant, while relative refinement makes the radial-distortion factors converge uniformly to one. Thus, the intrinsic criterion, together with relative refinement, provides the geometric controls used in anisotropic finite-element analysis without imposing a minimum-angle or shape-regularity condition.
Read the full abstract

Numerical Analysis

Heavy Ball GMRESR method for Nonsymmetric Linear Systems

The heavy ball GMRES (HBGMRES) method is one of Krylov subspace methods for linear systems combined with the restarted GMRES and the heavy ball method which is applied in optimization. HBGMRES not only keeps benefit of the restarted GMRES in limiting memory usage and controlling orthogonalization cost, but also is able to cover up the slow convergence problem in the restarted GMRES. Another type of Krylov subspace methods is GMRESR which consists of GCR as the outer algorithm and certain steps of GMRES as an inner method. Compared with HBGMRES, GMRESR gives the approximately optimal solution over the new search vector gained from GMRES and all previously kept search vectors in GCR. Even though GMRESR performs better than HBGMRES, it has slow convergence which is similar to GMRES and the restarted GMRES. Inspired by HBGMRES, we present the heavy ball GMRESR method (HBGMRESR) by using HBGMRES as the inner loop instead of using GMRES to salvage the lost convergence speed while still keeping the benefit of GMRESR that in the sense of a global minimization over some specific part of the Krylov subspace is done. Numerical tests on real data are presented to demonstratee the superiority of the new methods over GMRESR and HBGMRES.
Read the full abstract

Numerical Analysis

How different is rational approximation from piecewise polynomial approximation?

The first aim of this paper is to show that there is merit to the question posed in the title. Indeed, for certain function classes, approximation by rational functions and by piecewise polynomials are surprisingly similar. Between these two, polynomials are more widely studied and more widely used. One setting in which both approaches are equally well understood, and equally practical in use, is that of approximating univariate functions with point singularities. In this context we can fully address the question. We review classical literature on the topic which shows that both approaches do indeed achieve similar convergence rates. However, rational approximations come with significantly smaller constants. Owing to recent advances in practical rational approximation, we can augment the discussion with a comparison of modern numerical techniques that achieve the optimal rates. We end by showing numerically that the difference becomes even more pronounced in several variables.
Read the full abstract

Numerical Analysis

Randomized Tucker-Sketched GMRES

We address the problem of solving large-scale tensor-structured linear systems in the Tucker format. In this setting, standard iterative solvers such as GMRES face a fundamental bottleneck: the multilinear ranks of the Krylov basis vectors grow with the iteration count, leading to rapidly increasing tensor operation costs and memory requirements. To overcome these challenges, we propose two randomized algorithms within the sketched GMRES framework that replace full Arnoldi orthogonalization with short recurrences. The first, RHOSVD-Tucker sGMRES, uses randomized HOSVD with per-iteration rank selection, providing robustness across a wide range of problems. The second method, MLN-Tucker sGMRES, leverages the multilinear Nyström approximation with a fixed rank, enabling streaming computations; the streamability of the approximation further allows, at no additional cost, a memory-efficient reconstruction of the solution from a compact sketched representation of the Krylov basis. Both methods outperform standard low-rank Tucker solvers in symmetric and non-symmetric settings. Applied to inverse problems, the low-rank Tucker constraint acts as an implicit regularizer; combined with adaptive projected Tikhonov penalization and automatic regularization parameter selection, the methods yield stable reconstructions.
Read the full abstract

All authors Alberto Bucci, Martina Iannacito, Mirjeta Pasha, Rudi Smith

Numerical Analysis

Diffusion Quasi-Monte Carlo

We study high-dimensional numerical integration with respect to complex target measures using diffusion-based transport maps and randomized quasi-Monte Carlo (RQMC). Score-based diffusion models induce a deterministic probability flow ODE that transports a simple prior to the target, suggesting a principled way to transform low-discrepancy points on the unit cube into informative samples. We construct a cube-to-target map by composing a Gaussian base transformation (the component-wise inverse Gaussian CDF) with an Euler-discretized probability flow ODE. To retain unbiasedness under transport approximation, we formulate integration as importance sampling (IS) on the cube. Our main result provides verifiable conditions under which the resulting IS integrand satisfies the boundary growth condition, implying an $O(N^{-1+\epsilon})$ RMSE for scrambled nets. We then establish these conditions for diffusion probability-flow transport under mild bounded-derivative assumptions on the learned vector field, explicitly controlling the boundary singularities introduced by the inverse Gaussian CDF. Experiments range from a 2D mixture to 784D images and a 40,960D conditional vorticity-assimilation task; in the latter, blocked scrambled Sobol' sampling reduces the randomization standard deviation of nonlinear accuracy metrics at essentially unchanged online denoising cost. Together, these results give a theoretical and empirical foundation for combining diffusion generative modeling with high-precision RQMC integration.
Read the full abstract

Numerical Analysis · Cross-listed from cs.LG

Derivative Computation in PINNs: Automatic Differentiation, Finite Differences and Beyond

We systematically investigate finite-difference (FD) derivative computation in Physics-Informed Neural Networks (PINNs) as an alternative to automatic differentiation (AD). On three benchmark PDEs we show that, with a properly calibrated step size, FD matches AD in accuracy on every problem while running faster across the full tested batch-size range and using substantially less GPU memory, and that a stochastic variant we propose outperforms AD on a stationary problem. We further show that for neural architectures with inter-sample dependencies (e.g. BatchNorm, self-attention) the standard PyTorch autograd idiom is silently incorrect; the correct per-sample alternative is computationally infeasible at PINN-relevant batch sizes, while FD provides a forward-only approximation that is empirically an order of magnitude closer to the true per-sample derivative.
Read the full abstract

Numerical Analysis · Cross-listed from physics.geo-ph

A Spectral-Domain Pseudo-Inverse Method for True 3D Gravity Inversion

The main difficulty in 3D gravity inversion is that surface observations lack vertical wavenumber information, making the problem underdetermined and depth resolution poor. Building on the author's spectral-domain pseudo-inverse theory for unitary diagonalizable systems, this paper presents a true 3D inversion method. Surface data are analytically continued upward via Laplace's equation to form a 3D data volume, providing the vertical wavenumber sampling for the 3D Fourier transform. A general analytical expression for the half-space spectrum is derived, separating the horizontal spectrum from the vertical propagation kernel. The Green's function of the 3D Poisson equation is unitarily diagonalized, yielding the forward spectral response $\lambda(\kk) = -i 4\pi G k_z / K^2$, and a stable inverse filter $q_\alpha(\lambda) = \bar{\lambda}/(|\lambda|^2 + \alpha)$ is constructed. This filter is proved to have bounded stability and consistency (reducing to the exact inverse as $\alpha \to 0^+$). Validation with a homogeneous sphere model shows correct recovery of the anomaly location and singular behavior at the source. The contributions are twofold: (1) upward continuation constructs a 3D volume from 2D surface data, which does not add physical observations but unfolds the 3D spatial variation information (data innovation) derived through the governing physical equation; (2) the inversion is reduced to one forward 3D Fourier transform, one spectral scaling, and one inverse transform. The framework extends to any linear potential-field inverse problem with a translation-invariant forward operator, including magnetic inversion.
Read the full abstract

Numerical Analysis

Adaptive Time-Stepping Euler–Maruyama Scheme for SDEs with Non-Globally Lipschitz Coefficients: Uniform Convergence, Stability and Ergodicity

This paper develops an adaptive time-stepping Euler–Maruyama scheme for stochastic differential equations (SDEs) with non-globally Lipschitz drift and diffusion coefficients. By dynamically adjusting the timestep at each iteration, the proposed scheme effectively prevents numerical instability. We prove the moment boundedness of the numerical solution and establish a $1/2$-order strong convergence rate both on finite-time intervals and uniformly in time. Furthermore, the scheme faithfully inherits the $p$th moment exponential stability of the underlying SDE. For long-time ergodic dynamics, we establish the polynomial ergodicity of the numerical invariant measure. Moreover, we show that the numerical invariant measure converges to the invariant measure of the underlying SDE at an optimal rate of $1/2$ in the $L^q$-Wasserstein distance. Numerical experiments confirm our theoretical results and indicate the superior accuracy and computational performance of the proposed scheme over several fixed-step and existing adaptive methods.
Read the full abstract

Numerical Analysis

Locking analysis and high-order locking-free finite element method for three-dimensional electroporoelasticity equations

Electroporoelasticity equations couple Maxwell's equations with Biot's poroelasticity model and admit severe Poisson locking in standard conforming finite element discretizations when the Lamé constant is large. In this paper, we provide a rigorous analysis of the Poisson locking phenomenon for three-dimensional quasi-static electroporoelasticity equations and show that the spatial convergence order of conforming finite element approximations is reduced in the nearly incompressible regime. To eliminate this locking effect, we introduce a five-field formulation and a fully discrete high-order finite element method. We prove uniform stability with respect to large Lamé constant. Based on this, we establish a locking-free scheme by deriving its uniform error estimates with respect to the Lamé coefficient. The analysis covers the fully coupled electromagnetic-poroelastic system and applies to high-order elements in three dimensions. Extensive numerical experiments are presented to verify the theoretical convergence rates and demonstrate robustness with respect to the Lamé parameter.
Read the full abstract

Numerical Analysis

High-order $L^2$-Galerkin schemes for the decoupled potential integral equations

We prove Sobolev space well-posedness of the decoupled potential integral equations. Using basic pseudo-differential techniques, discrete stability is obtained for the standard $L^2$-pairing, and we deduce spectral convergence for a high-order scheme employing discontinuous basis functions, or with additional element-wise tangential continuity imposed in the vectorial part of the system. The scheme is shown to be much better conditioned than a mature commercial EFIE/CFIE solver. In many cases appearing almost independent of wavenumber, mesh density, and the order.
Read the full abstract

Numerical Analysis

Successive Schur-Riesz Analysis for Approximation

Many approximation methods enlarge a trial space by adjoining function blocks generated by different operators. Exact redundancy and strong cross-level interaction can make coefficients nonunique and render pairwise or diagonal-dominance tests needlessly pessimistic. For \(V_m=\sum_{\ell\leq m}S_\ell(E_\ell)\) in a Hilbert space $\mathcal H$, we quotient coefficients representing the same function and control successive orthogonal innovations to obtain Riesz bounds independent of \(m\). The setting includes factored operators \(S_\ell=T_\ell\circ\cdots\circ T_1:E_\ell\to\mathcal H\), with compatible intermediate spaces, but the theorem allows arbitrary bounded \(S_\ell\). The same constants control approximation, truncation, perturbation, and levelwise error. A block Schur complement identifies the intrinsic new dimension and gives the exact reduction in squared best-approximation error, leading to a constructive enrichment procedure. Nonstationary and lifted examples give positive intrinsic bounds where diagonal-dominance estimates are negative or labelled Gram matrices are singular; adaptive and recycled-subspace calculations illustrate the distinct roles of representation stability and application-specific utility.
Read the full abstract

Numerical Analysis · Cross-listed from cs.LG

Long-Time Trajectory Approximation via SA-NODEs: Model Predictive and Floquet Strategies

We study the approximation of dynamical systems by semi-autonomous neural ordinary differential equations (SA-NODEs) over long time horizons. For a single network trained on the whole horizon, the available error bound deteriorates double exponentially in the horizon length. We develop two training strategies that avoid this barrier, each built on a reset of the state. The model predictive strategy partitions the horizon adaptively and restarts every window from observed data: when training meets a prescribed tolerance on every window, the composite model meets it uniformly in time, with a parameter budget linear in the horizon for targets with a bounded, uniformly regular reachable tube. The Floquet strategy addresses autonomous targets with a stable limit cycle and uses no data at deployment: a certified contraction of the learned return map confines the error to linear growth in the number of elapsed periods. For the time-periodic architecture we deploy, the scalar certificate degenerates; we prove instead a uniform-in-time orbital guarantee whose hypotheses are measured on the trained model, and an obstruction showing that, for an exactly periodic learned field, small one-period error and a contracting stroboscopic map cannot hold at once. Numerical experiments on four benchmarks confirm the predicted error laws and measure the hypotheses of every guarantee.
Read the full abstract