Society for Industrial and Applied Mathematics: SIAM Journal on Scientific Computing: Table of Contents

A Null Infinity Layer for Wave Scattering

Anil Zenginoğlu — 2026-05-04T07:00:00Z

SIAM Journal on Scientific Computing, Volume 48, Issue 3, Page A1075-A1100, June 2026.
Abstract. We solve time-harmonic wave scattering problems on unbounded domains without domain truncation by mapping the unbounded domain to a bounded domain and scaling the oscillatory decay towards infinity. The technique, first developed in numerical relativity for time-domain wave equations, solves for the far-field pattern using compactification at infinity, avoiding the outer boundary problem. We design a layer that restricts the transformations to an annular domain. The resulting null infinity layer solves both the outer boundary and radiation extraction problems. We demonstrate its practical implementation using finite difference, Chebyshev collocation, spectral Galerkin, and finite element methods in one and two dimensions.

Generalized Moving Least-Squares Methods for Solving Vector-Valued PDEs on Unknown Manifolds

Rongji Li — 2026-05-04T07:00:00Z

SIAM Journal on Scientific Computing, Volume 48, Issue 3, Page A1101-A1132, June 2026.
Abstract. In this paper, we extend the generalized moving least-squares method in two different ways to solve the vector-valued PDEs on unknown smooth manifolds without boundaries, identified with randomly sampled point cloud data. The two approaches are referred to as the intrinsic method and the extrinsic method. For the intrinsic method, which relies on local approximations of metric tensors, we simplify the formula of Laplacians and covariant derivatives acting on vector fields at the base point by calculating them in a local Monge coordinate system. On the other hand, the extrinsic method formulates tangential derivatives on a submanifold as the projection of the directional derivative in the ambient Euclidean space onto the tangent space of the submanifold. One challenge of this method is that the discretization of vector Laplacians yields a matrix whose size relies on the ambient dimension. To overcome this issue, we reduce the dimension of vector Laplacian matrices by employing a coordinate transformation. The complexity of both methods scales well with the dimension of manifolds rather than with the ambient dimension. We also present supporting numerical examples, including eigenvalue problems, linear Poisson equations, and nonlinear Burgers’ equations, to examine the numerical accuracy of the proposed methods on various smooth manifolds.

Solving Eigenvalue Problems in High Dimensions Using Contour Integration and Tensor Train Format

Olivier Coulaud — 2026-05-04T07:00:00Z

SIAM Journal on Scientific Computing, Volume 48, Issue 3, Page A1133-A1152, June 2026.
Abstract. In high-dimensional settings, solving eigenvalue problems is hindered by the curse of dimensionality, particularly when only a subset of eigenpairs within a prescribed spectral interval is sought. In this work, we investigate an adaptation of the FEAST algorithm, originally developed for symmetric eigenproblems based on contour integration, to computations where both operators and vectors are represented in the tensor train (TT) format. This representation drastically reduces memory and computational demands. We introduce an adaptive scheme for determining the projection subspace dimension by incorporating a rank-revealing Modified Gram–Schmidt procedure with pivoting tailored to TT-vectors. A perturbation-based analysis provides explicit bounds on the attainable residual accuracy, from which we derive a robust stopping criterion for the proposed TT-FEAST algorithm. Moreover, we design a continuation strategy that gradually refines convergence and rounding tolerances to effectively control memory growth during iterations. To demonstrate the effectiveness of TT-FEAST as a viable alternative to existing high-dimensional eigensolvers when a few eigenvalues are required, we present numerical experiments on problems up to 12 dimensions, including the Laplacian and a vibrational Hamiltonian operator.

On a Shrink-and-Expand Technique for Symmetric Block Eigensolvers

Yuqi Liu — 2026-05-04T07:00:00Z

SIAM Journal on Scientific Computing, Volume 48, Issue 3, Page A1153-A1183, June 2026.
Abstract. In block eigenvalue algorithms, such as the subspace iteration algorithm and the locally optimal block preconditioned conjugate gradient (LOBPCG) algorithm, a large block size is often employed to achieve robustness and rapid convergence. However, using a large block size also increases the computational cost. Traditionally, the block size is typically reduced after convergence of some eigenpairs, known as deflation. In this work, we propose a nondeflation-based, more aggressive technique, where the block size is adjusted dynamically during the algorithm. This technique can be applied to a wide range of block eigensolvers, reducing computational cost without compromising convergence speed. We present three adaptive strategies for adjusting the block size, and apply them to four well-known eigensolvers as examples. Theoretical analysis and numerical experiments are provided to illustrate the efficiency of the proposed technique. In practice, an overall acceleration of 20% to 30% is observed. Reproducibility of computational results. This paper has been awarded the “SIAM Reproducibility Badge: Code and data available” as a recognition that the authors have followed reproducibility principles valued by SISC and the scientific computing community. Code and data that allow readers to reproduce the results in this paper are available at https://github.com/IzayoiYuuki/Shrink-and-expand_public and in the supplementary materials (Shrink-and-expand_public-master.zip [70.2MB]).

Solving Random Hyperbolic Conservation Laws Using Linear Programming

Shaoshuai Chu — 2026-05-05T07:00:00Z

SIAM Journal on Scientific Computing, Volume 48, Issue 3, Page A1184-A1205, June 2026.
Abstract. A novel structure-preserving numerical method to solve random hyperbolic systems of conservation laws is presented. The method uses a concept of generalized, measure-valued solutions to random conservation laws. This yields a linear partial differential equation with respect to the Young measure and allows for the computation of the approximation based on linear programming problems. We analyze structure-preserving properties of the derived numerical method and discuss its advantages and disadvantages. We numerically demonstrate the approach on the one-dimensional Burgers and isentropic Euler equations and compare with stochastic collocation. In addition, we introduce a discontinuous-flux test in which different entropies used in the linear-program objective select different weak entropy solutions, and we report the corresponding changes in the moments and supports of the Young measure.

A Variable Dimension Sketching Strategy for Nonlinear Least-Squares

Stefania Bellavia — 2026-05-05T07:00:00Z

SIAM Journal on Scientific Computing, Volume 48, Issue 3, Page A1206-A1234, June 2026.
Abstract. We present a stochastic inexact Gauss–Newton method for the solution of nonlinear least-squares. To reduce the computational cost with respect to the classical method, at each iteration the proposed algorithm approximately minimizes the local model on a random subspace. The dimension of the subspace varies along the iterations, and two strategies are considered for its update: the first is based solely on the Armijo condition, and the latter is based on information from the true Gauss–Newton model. Under suitable assumptions on the objective function and the random subspace, we prove a probabilistic bound on the number of iterations needed to drive the norm of the gradient below any given threshold. Moreover, we provide a theoretical analysis of the local behavior of the method. The numerical experiments demonstrate the effectiveness of the proposed method.

An Energy-Stable Parametric Finite Element Method for Willmore Flow with Normal-Tangential Velocity Splitting

Harald Garcke — 2026-05-05T07:00:00Z

SIAM Journal on Scientific Computing, Volume 48, Issue 3, Page A1235-A1259, June 2026.
Abstract. We propose and analyze an energy-stable fully discrete parametric approximation for Willmore flow of hypersurfaces in two and three space dimensions. We allow for the presence of spontaneous curvature effects and for open surfaces with boundary. The presented scheme is based on a new geometric partial differential equation that combines an evolution equation for the mean curvature with a separate equation that prescribes the tangential velocity. The mean curvature is used to determine the normal velocity within the gradient flow structure, thus guaranteeing unconditional energy stability for the discrete solution upon suitable discretization. We introduce a novel weak formulation for this geometric PDE, in which different types of boundary conditions can be naturally enforced. We further discretize the weak formulation to obtain a fully discrete parametric finite element method, for which well-posedness can be rigorously shown. Moreover, the constructed scheme admits an unconditional stability estimate in terms of the discrete energy. Extensive numerical experiments are reported to showcase the accuracy and robustness of the proposed method for computing Willmore flow of both curves in [math] and surfaces in [math].

Evaluation of Resonances: Adaptivity and AAA Rational Approximation of Randomly Scalarized Boundary Integral Resolvents

Oscar P. Bruno — 2026-05-06T07:00:00Z

SIAM Journal on Scientific Computing, Volume 48, Issue 3, Page A1260-A1283, June 2026.
Abstract. This paper introduces a novel algorithm that, employing rational approximants of randomly scalarized boundary integral resolvents, efficiently evaluates acoustic and electromagnetic resonances in both open and closed cavities. The desired cavity resonances (also known as “eigenvalues” for interior problems and “scattering poles” or “complex eigenvalues” for exterior and open-cavity problems) are obtained as the poles of associated rational approximants; both the approximants and their poles are produced by means of the recently introduced AAA rational approximation algorithm. In fact, the proposed resonance search method applies to any nonlinear eigenvalue problem associated with a given function [math], wherein, denoting [math], a complex value [math] is sought for which [math] for some nonzero [math]. For the scattering problems considered in this paper, which include interior, exterior, and open-cavity problems, [math] is taken to equal a spectrally discretized version of a Green function–based boundary integral operator at spatial frequency [math]. In all cases, the scalarized resolvent is given by an expression of the form [math], where [math] are fixed random vectors. The proposed adaptive search strategy relies on use of a rectangular subdivision of the resonance search domain which is locally refined to ensure that all resonances in the domain are captured. The approach works equally well in the case in which the search domain is a one-dimensional set, such as, e.g., an interval of the real line, in which case the rectangles used degenerate into subintervals of the search domain. A variety of numerical results are presented, including comparisons with well-known methods based on complex contour integration, and a discussion of the asymptotics that result as open cavities approach closed cavities—in all, demonstrating the accuracy provided by the method, for low- and high-frequency states alike.

Energy Dissipating ALE-MDR Method for Navier–Stokes Free Boundary Problems with Moving Contact Line

Jiashun Hu — 2026-05-06T07:00:00Z

SIAM Journal on Scientific Computing, Volume 48, Issue 3, Page A1284-A1311, June 2026.
Abstract. We propose a novel arbitrary Lagrangian–Eulerian (ALE) finite element method for simulating incompressible Navier–Stokes flows driven by surface tension on evolving free boundaries as well as for simulating the dynamics of moving contact lines at fluid–solid interfaces under a prescribed contact angle condition. The proposed method ensures energy dissipation while incorporating the following two key strategies to maintain mesh quality within the ALE framework: the artificial tangential motion strategy of Barrett, Garcke, and Nürnberg (BGN) on the fluid’s free surface, and the minimal-deformation-rate (MDR) approach within the fluid bulk and along the solid boundary. A central challenge in three-dimensional simulations—the ambiguity of tangential mesh velocity at the moving contact line—is resolved by enforcing a constraint that eliminates tangential motion at the contact line, thereby maintaining stability of tangential motion and mesh quality of computed surfaces. Numerical experiments in both two and three dimensions demonstrate the robustness of the proposed method in maintaining mesh quality and energy dissipation.

An Additive Two-Level Parallel Variant of the DMRG Algorithm with Coarse-Space Correction

Laura Grigori — 2026-05-07T07:00:00Z

SIAM Journal on Scientific Computing, Volume 48, Issue 3, Page A1312-A1337, June 2026.
Abstract. The density matrix renormalization group (DMRG) algorithm is a popular alternating minimization scheme for solving high-dimensional optimization problems in the tensor train format. Classical DMRG, however, is based on sequential minimization, which raises challenges in its implementation on parallel computing architectures. To overcome this, we propose a novel additive two-level DMRG algorithm that combines independent, local minimization steps with a global update step using a subsequent coarse-space minimization. Our proposed algorithm, which is directly inspired by additive Schwarz methods from the domain decomposition literature, is particularly amenable to implementation on parallel, distributed architectures since both the local minimization steps and the construction of the coarse-space can be performed in parallel. Numerical experiments on strongly correlated molecular systems demonstrate that the method achieves competitive convergence rates while achieving significant parallel speedups.

A Structure-Preserving, Helicity-Conserving High-Order ALE Finite Element Method for Compressible Ideal MHD Systems

Shipeng Mao — 2026-05-08T07:00:00Z

SIAM Journal on Scientific Computing, Volume 48, Issue 3, Page A1338-A1366, June 2026.
Abstract. Structure-preservation, accuracy, and robustness are critical aims for magnetohydrodynamic (MHD) simulations. In this paper, we propose a high-order, magnetic helicity conserving, divergence-free, and positivity-preserving arbitrary-Lagrangian-Eulerian (ALE) finite element method for compressible ideal MHD systems. Notably, helicity conservation can be achieved using arbitrary high-order time integration methods, overcoming the limitations of Crank–Nicolson methods commonly used in existing literature. To our knowledge, this is the first Lagrangian and ALE scheme that addresses magnetic helicity conservation in MHD. These structure-preserving properties arise from a new Lagrangian MHD formulation with differential forms and an extended definition of material derivatives through Lie derivatives. By utilizing high-order finite element exterior calculus (FEEC) discretization, these properties are maintained at negligible cost during Lagrangian steps. To overcome mesh distortion issues, we introduce a mesh optimization strategy that automatically smooths the mesh based on quality indicators. After mesh smoothing steps, we propose structure-preserving remapping methods, including a positivity-preserving remap for density and a high-order vector potential-based remap for magnetic fields that preserves both divergence-free condition and magnetic helicity.

An [math] Helmholtz Solver Using WaveHoltz and Overset Grids

Daniel Appelö — 2026-05-14T07:00:00Z

SIAM Journal on Scientific Computing, Volume 48, Issue 3, Page A1367-A1398, June 2026.
Abstract. We develop efficient and high-order accurate solvers for the Helmholtz equation on complex geometry. The solvers are based on the WaveHoltz algorithm, which computes solutions of the Helmholtz equation by time-filtering solutions of the wave equation. The approach avoids the need to invert an indefinite matrix which can cause convergence difficulties for many iterative solvers for indefinite Helmholtz problems. Complex geometry is treated with overset grids which use Cartesian grids throughout most of the domain together with curvilinear grids near boundaries. The basic WaveHoltz fixed-point iteration is accelerated using GMRES and also by a deflation technique using a set of precomputed eigenmodes. The solution of the wave equation is solved efficiently with implicit time-stepping using as few as five time-steps per period, independent of the mesh size. When multigrid is used to solve the implicit time-stepping equations, the cost of the resulting WaveHoltz scheme scales linearly with the total number of grid points [math] (at fixed frequency) and is thus optimal in CPU time and memory usage as the mesh is refined. Numerical results are given for problems in two and three dimensions, to second- and fourth-order accuracy, and they show the potential of the approach to solve a wide range of large-scale problems.

A Stable Matrix Version of the Wideband Fast Multipole Method for the 2D Helmholtz Kernel

Michelle Michelle — 2026-05-27T07:00:00Z

SIAM Journal on Scientific Computing, Volume 48, Issue 3, Page A1399-A1429, June 2026.
Abstract. This paper presents a stable matrix version of the wideband fast multipole method (FMM) for the 2D Helmholtz kernel. It is known that the FMM may experience stability issues in both high-frequency and low-frequency regimes, some of which can be mitigated and others are inherent in nature. Inspired by recent studies, we propose a balancing strategy to overcome the stability challenge that exists in the low-frequency regime. The balancing strategy utilizes some simple properties of Bessel and Hankel functions so as to produce theoretically guaranteed norm bounds for relevant low-rank expansion factors and translation operators. We then present an intuitive and stable matrix version of the wideband FMM, which utilizes two different expansions of the 2D Helmholtz kernels: one that always behaves well in the low-frequency regime based on our balancing strategy, and the other that behaves well (under certain conditions) in the high-frequency regime. The backward stability of this wideband FMM is rigorously justified based on our studies of the norm bounds of the low-rank factors and translation operators. Some numerical experiments demonstrate the effectiveness and the accuracy of the wideband FMM.

Fast and Accurate Computation of Classical Gaussian Quadratures

Amparo Gil — 2026-05-07T07:00:00Z

SIAM Journal on Scientific Computing, Volume 48, Issue 3, Page B289-B316, June 2026.
Abstract. Algorithms for computing the classical Gaussian quadrature rules (Gauss–Jacobi, Gauss–Laguerre, and Gauss–Hermite) are presented, based on globally convergent fourth order iterative methods combined with asymptotic approximations, which are applied in complementary regions of the parameter space. This approach yields methods that improve upon existing algorithms in speed, accuracy, and computational range. The MATLAB algorithm for Gauss–Jacobi is faster than previous methods and lifts the upper restrictions on the parameters imposed by those methods ([math]); for example, for degrees up to [math] all nodes and weights can be computed within the underflow limit for [math], and the computable range of parameters is much larger for smaller degrees, limited only by intrinsic overflow/underflow constraints. For the particular case of Gauss–Legendre quadrature ([math]), a specific asymptotic approach is considered, which yields the most efficient MATLAB implementation available so far. The Gauss–Laguerre and Gauss–Hermite algorithms incorporate subsampling, and scaling is also available in order to extend the computational range. Gauss–Radau and Gauss–Lobatto variants are also considered, along with the computation of the associated barycentric weights. Additionally, arbitrary-precision algorithms (in Maple) are offered for the symmetric cases (Gauss–Gegenbauer and Gauss–Hermite), which can be used to compute thousands of nodes with hundreds of digits in a matter of seconds.

Parallel Distributed Out-of-Core Solvers with Low-Rank Compression for Coupled FEM/BEM Systems

Emmanuel Agullo — 2026-05-08T07:00:00Z

SIAM Journal on Scientific Computing, Volume 48, Issue 3, Page B317-B336, June 2026.
Abstract. Distributed-memory parallelism allows sparse and dense direct solvers to process ever larger linear systems. While some of them also leverage out-of-core computation and low-rank compression to reduce memory footprint and computation time, the use of these techniques for solving coupled finite elements method FEM/BEM linear systems that contain both sparse and dense parts has been little investigated. In [], we therefore introduced two classes of algorithms, namely, the multi-solve and multi-factorization algorithms, which use low-rank compression to solve relatively large coupled FEM/BEM systems on a shared-memory machine. Nevertheless, to the best of our knowledge, no existing approaches exploit out-of-core computation, combine it with low-rank compression, or apply these techniques in distributed memory. In this paper, we propose a design of the multi-solve and multi-factorization algorithms that, beyond low-rank compression, incorporate out-of-core computation and distributed-memory parallelism to process very large coupled FEM/BEM systems. An experimental study on up to 16 computation nodes, each equipped with 48 cores and 180 GB of RAM, shows that the proposed algorithms, implemented on top of state-of-the-art sparse and dense direct solvers, can process systems of up to 42 million and 7 million unknowns, respectively, instead of 6 million unknowns with a standard sparse/dense solver coupling.

Unified Interface Flux Evaluation in a General Discontinuous Galerkin Spectral Element Framework

Boyang Xia — 2026-05-14T07:00:00Z

SIAM Journal on Scientific Computing, Volume 48, Issue 3, Page B337-B359, June 2026.
Abstract. High-order discontinuous Galerkin spectral element methods (DGSEM) have received growing attention and development in recent years, especially in the regime of computational fluid dynamics. The inherent flexibility of the DG approach in handling nonconforming interfaces, such as those encountered in moving geometries or [math]-refinement, presents a significant advantage for real-world simulations. Despite the well-established mathematical framework of DG methods, practical implementation challenges persist in boosting performance and capability. Most previous studies only focus on certain choices of element shape or basis type in a structured mesh, although they have demonstrated the capability of DGSEM in complex flow simulations. This work discusses low-cost and unified interface flux evaluation approaches for general spectral elements in unstructured meshes, alongside their implementations in the open-source spectral element framework, Nektar++. The initial motivation arises from the discretization of Helmholtz equations by the symmetric interior penalty method, in which the system matrix can easily become nonsymmetric if the flux is not properly evaluated on nonconforming interfaces. We focus on the polynomial nonconforming case in this work, but extending to the geometric nonconforming case is theoretically possible. Comparisons of different approaches, tradeoffs, and performance of our initial matrix-free implementation are also included, contributing to the broader discourse on high-performance spectral element method implementations.

A Massively Parallel Interior-Point Method for Arrowhead Linear Programs with Local Linking Structure

Nils–Christian Kempke — 2026-05-21T07:00:00Z

SIAM Journal on Scientific Computing, Volume 48, Issue 3, Page B360-B385, June 2026.
Abstract. In practice, nonspecialized interior-point algorithms often cannot utilize the massively parallel compute resources offered by modern many- and multicore compute platforms. However, efficient distributed solution techniques are required, especially for large-scale linear programs. This article describes a new decomposition technique for systems of linear equations implemented in the parallel interior-point solver PIPS-IPM++. The algorithm exploits a matrix structure commonly found in optimization problems: a doubly bordered block-diagonal or arrowhead structure with linking constraints and variables often only linking few, consecutive blocks. This structure is preserved in the linear Karush-Kuhn-Tucker (KKT) systems solved during each iteration of the interior-point method. We present a hierarchical Schur complement decomposition that distributes and solves the linear optimization problem. It is designed for high-performance architectures and scales well with the availability of additional computing resources. The decomposition approach uses the border constraints’ locality to decouple the factorization process. Our approach is motivated by large-scale economic dispatch problems but can also be applied to other problem classes. We demonstrate the performance of our method on a set of mid- to large-scale instances, some of which have more than [math] nonzeros in their constraint matrices.

An Implementation of Quantum Oracles for the Finite Element Method

Sven Danz — 2026-06-01T07:00:00Z

SIAM Journal on Scientific Computing, Volume 48, Issue 3, Page B386-B419, June 2026.
Abstract. In order to assess potential advantages of quantum algorithms that require quantum oracles as subroutines, the careful evaluation of the overall complexity of the oracles themselves is crucial. This study examines the quantum routines required for the implementation of oracles used in the block-encoding of the [math] stiffness and mass matrices, which typically emerge in the finite element analysis of elastic structures. Starting from basic quantum adders, we show how to construct the necessary oracles, which require the calculation of polynomials, square root, and the implementation of conditional operations. We propose quantum subroutines based on fixed-point arithmetic that, given an [math]-qubit register, construct the oracle using [math] ancilla qubits and have an [math] runtime, with [math] the order at which we truncate the polynomials and [math] the number of iterations in the Newton–Raphson subroutine for the square root, while [math] and [math] are the number of hypercuboids used to approximate the geometry and the boundary, respectively. Since in practice [math] scales as [math], for numbers between 0 and [math], and assuming that the other parameters are fixed independently of [math], this shows that the oracles, while still costly in practice, do not endanger potential polynomial or exponential advantages in [math].

LREI: A Fast Numerical Solver for Quantum Landau–Lifshitz Equations

Davoud Mirzaei — 2026-06-01T07:00:00Z

SIAM Journal on Scientific Computing, Volume 48, Issue 3, Page B420-B439, June 2026.
Abstract. We develop LREI, short for low-rank eigenmode integration, a memory- and time-efficient numerical scheme for solving quantum Landau–Lifshitz (q-LL) and quantum Landau–Lifshitz–Gilbert (q-LLG) equations, which govern spin dynamics in open quantum systems. Although the system size grows exponentially with the number of spins, our approach benefits from the low-rank structure of the density matrix and the sparsity of system Hamiltonians to avoid costly full matrix computations. By representing density matrices in terms of their low-rank factors and using Krylov subspace techniques for partial eigendecompositions, we reduce the per-step complexity of Runge–Kutta and Adams–Bashforth schemes from [math] to [math], where [math] is the Hilbert space dimension for [math] spins and [math] is the effective rank of the density matrix. Likewise, the memory footprint is reduced from [math] to [math], since no full [math] matrices are ever formed. Among several technical improvements we applied, one key idea was to handle the computation of the action of the invariant subspace of the density matrix associated with its zero eigenvalues. This was accomplished by applying Householder reflectors constructed for the dominant eigenspace, thereby enabling the entire solution process to proceed without ever forming any large matrices. As an example, we can now evolve a time step of a 20-spin system—corresponding to a density matrix size exceeding one million—in just a few seconds on a standard laptop. In addition, both classes of Runge–Kutta and Adams–Bashforth techniques are reformulated to preserve the physical properties of the density matrix throughout the time evolution. The new low-rank algorithm makes it possible to simulate much larger spin systems that were previously computationally infeasible. This, in turn, provides a powerful tool for comparing the q-LL and q-LLG dynamics, assessing the validity of each model, and exploring how quantum features such as correlations and entanglement evolve across different regimes of system size and damping.

Neural Networks for Bayesian Inverse Problems Governed by a Nonlinear ODE

German Villalobos — 2026-05-04T07:00:00Z

SIAM Journal on Scientific Computing, Volume 48, Issue 3, Page C415-C452, June 2026.
Abstract. We investigate the use of neural networks (NNs) for the estimation of hidden model parameters and uncertainty quantification from noisy observational data for inverse parameter estimation problems. We formulate the parameter estimation as a Bayesian inverse problem. We consider a parametrized system of nonlinear ordinary differential equations (ODEs), which is the FitzHugh–Nagumo model. The considered problem exhibits significant mathematical and computational challenges for classical parameter estimation methods, including strong nonlinearities, nonconvexity, and sharp gradients. We explore how NNs overcome these challenges by approximating reconstruction maps for parameter estimation from observational data. The considered data are time series of the spiking membrane potential of a biological neuron. We infer parameters controlling the dynamics of the model, noise parameters of autocorrelated additive noise, and noise modeled via stochastic differential equations, as well as the covariance matrix of the posterior distribution to expose parameter uncertainties—all with just one forward evaluation of an appropriate NN. We report results for different NN architectures and study the influence of noise on prediction accuracy. We also report timing results for training NNs on dedicated hardware. Our results demonstrate that NNs are a versatile tool to estimate parameters of the dynamical system, stochastic processes, as well as uncertainties as they propagate though the governing ODE.

Rank-Inspired Neural Network for Solving Partial Differential Equations

Wentao Peng — 2026-05-06T07:00:00Z

SIAM Journal on Scientific Computing, Volume 48, Issue 3, Page C453-C478, June 2026.
Abstract. This paper proposes a rank-inspired neural network (RINN) to tackle the initialization sensitivity issue of physics-informed extreme learning machines (PIELMs) when numerically solving PDEs. Unlike PIELM, which randomly initializes the parameters of its hidden layers, RINN incorporates a preconditioning stage. In this stage, covariance-driven regularization is employed to optimize the orthogonality of the basis functions generated by the last hidden layer. The key innovation lies in minimizing the off-diagonal elements of the covariance matrix derived from the hidden-layer output. By doing so, pairwise orthogonality constraints across collocation points are enforced, which effectively enhances both the numerical stability and the approximation ability of the optimized function space. The RINN algorithm unfolds in two sequential stages. First, it conducts a nonlinear optimization process to orthogonalize the basis functions. Subsequently, it solves the PDE constraints using the linear least-squares method. Extensive numerical experiments demonstrate that RINN significantly reduces performance variability due to parameter initialization compared to PIELM. Incorporating an early stopping mechanism based on PDE loss further improves stability, ensuring consistently high accuracy across diverse initialization settings. Reproducibility of computational results. This paper has been awarded the “SIAM Reproducibility Badge: Code and data available” as a recognition that the authors have followed reproducibility principles valued by SISC and the scientific computing community. Code and data that allow readers to reproduce the results in this paper are available at https://github.com/Lambda-None/RINN and in the supplementary materials (RINN-main.zip [5.63MB]).

PODNO: Proper Orthogonal Decomposition Neural Operators

Zilan Cheng — 2026-05-11T07:00:00Z

SIAM Journal on Scientific Computing, Volume 48, Issue 3, Page C479-C504, June 2026.
Abstract. In this paper, we introduce proper orthogonal decomposition neural operators (PODNO) for solving partial differential equations (PDEs) dominated by high-frequency components. Building on the structure of Fourier neural operators (FNO), PODNO replaces the Fourier transform with (inverse) orthonormal transforms derived from the proper orthogonal decomposition (POD) method to construct the integral kernel. Due to the optimality of POD basis, the PODNO has the potential to outperform FNO in both accuracy and computational efficiency for high-frequency problems. From an analysis point of view, we established the universality of a generalization of PODNO, termed generalized spectral operators. In addition, we evaluate PODNO’s performance numerically on dispersive equations such as the nonlinear Schrödinger equation and the Kadomtsev–Petviashvili equation. Reproducibility of computational results. This paper has been awarded the “SIAM Reproducibility Badge: Code and data available” as a recognition that the authors have followed reproducibility principles valued by SISC and the scientific computing community. Code and data that allow readers to reproduce the results in this paper are available at https://github.com/Zilan-Cheng/PODNO and in the supplementary material (PODNO-main.zip PODNO-main.zip [45 KB]).

Machine Learning-Based Quadratic Closures for Non-Intrusive Reduced Order Models

Gabriele Codega — 2026-05-14T07:00:00Z

SIAM Journal on Scientific Computing, Volume 48, Issue 3, Page C505-C525, June 2026.
Abstract. In the present work, we introduce a data-driven approach to enhance the accuracy of non-intrusive reduced order models (ROMs). In particular, we focus on ROMs built using proper orthogonal decomposition (POD) in an under-resolved and marginally-resolved regime, i.e., when the number of modes employed is not enough to capture the system dynamics. We propose a method to re-introduce the contribution of neglected modes through a quadratic correction term, given by the action of a quadratic operator on the POD coefficients. Differently from the state-of-the-art methodologies, where the operator is learned via least-squares optimization [R. Geelen, S. Wright, and K. Willcox, Comput. Methods Appl. Mech. Engrg., 403 (2023), 115717; J. Barnett and C. Farhat, J. Comput. Phys., 464 (2022), 111348], we propose to parametrize the operator by a Multi-Input Operators Network (MIONet). This way, we are able to build models with higher generalization capabilities, where the operator itself is continuous in space—thus agnostic of the domain discretization—and parameter-dependent. We test our model on two standard benchmarks in fluid dynamics and show that the correction term improves the accuracy of standard POD-based ROMs.

A Direct Method for Computing the Complex-Valued Triple Decomposition of Third-Order Tensors

Yannan Chen — 2026-03-04T08:00:00Z

SIAM Journal on Scientific Computing, Volume 48, Issue 2, Page A425-A447, April 2026.
Abstract. A direct method is first proposed to compute a complex-valued triple decomposition of a third-order tensor. The method operates under three assumptions: (i) The triple rank [math] of an [math] tensor satisfies [math]. (ii) Two factor tensors of the triple decomposition are generic. (iii) The third factor tensor has linearly independent fibers. When [math], the computational cost of the proposed direct method is approximately [math] flops. Moreover, a sufficient condition for the essential uniqueness of the tensor triple decomposition is established under these assumptions. Numerical experiments demonstrate that our direct method is at least 10 times faster than alternating least squares and optimization-based iterative approaches. Finally, we showcase applications of the method in large-scale video analysis and stochastic partial differential equations using triple tensor representations.

Computation of Shape Taylor Expansions for Scattering Problems in Two Dimensions

Gang Bao — 2026-03-04T08:00:00Z

SIAM Journal on Scientific Computing, Volume 48, Issue 2, Page A448-A473, April 2026.
Abstract. Shape derivative is an important analytical tool for studying scattering problems involving perturbations in scatterers. Many applications, including inverse scattering, optimal design, and uncertainty quantification, are based on shape derivatives. However, computing high order shape derivatives is challenging due to the complexity of shape calculus. This work introduces a comprehensive method for computing shape Taylor expansions in two dimensions using recurrence formulas. The approach is developed under sound-soft, sound-hard, impedance, and transmission boundary conditions. Additionally, we apply the shape Taylor expansion to uncertainty quantification in wave scattering, enabling high order moment estimation for the scattered field under random boundary perturbations. Numerical examples are provided to illustrate the effectiveness of the shape Taylor expansion in achieving high order approximations.

A Class of High Order Unconditional Maximum Principle Preserving Extended IFRK Schemes for Allen-Cahn Equation

Zhuhan Jiang — 2026-03-06T08:00:00Z

SIAM Journal on Scientific Computing, Volume 48, Issue 2, Page A474-A493, April 2026.
Abstract. This paper proposes a novel framework for constructing up to sixth-order, explicit, and unconditional maximum principle preserving extended integrating factor Runge–Kutta (eIFRK) schemes for the Allen-Cahn equation. The proposed framework significantly relaxes the constraints on the exponential coefficients and offers a recursive procedure to establish high-order schemes. Satisfying the requirement of the local truncation error in the numerical integration, derived from Duhamel’s formula, provides an innovative approach to achieve high-order accuracy instead of solving the order conditions directly. Specifically, in order to derive the target order scheme, Gauss quadrature is applied to approximate the implicit integration, where one order lower scheme is employed to compute explicitly at quadrature points. Consequently, maximum bound principle and [math] norm analysis are uniformly proved for eIFRK4–6. Numerical experiments on eIFRK schemes validate their accuracy, maximum-norm preservation, and energy stability. In comparison to the classic second- and third-order sIFRK methods [SIAM J. Sci. Comput., Vol 43, 2021, pp. A1780-A1802], the presented high-order eIFRK schemes, particularly eIFRK6, exhibit enhanced efficiency and stability during long-term simulations, especially for large time steps.

Differential Inversion of the Implicit Euler Method

Uwe Naumann — 2026-03-06T08:00:00Z

SIAM Journal on Scientific Computing, Volume 48, Issue 2, Page A494-A511, April 2026.
Abstract. The implicit Euler method integrates systems of ordinary differential equations [math] with differentiable right-hand side [math] from an initial state [math] to a target time [math] as [math]. Discretization of the time interval [math] yields [math] time steps. We present a method for efficiently computing the product of its inverse Jacobian [math] with a given vector [math]. We show that the differential inverse [math] can be evaluated for given [math] with a computational cost of [math] as opposed to the widely used [math] or, if applying Algorithmic Differentiation to [math] naively, even [math]. The theoretical results are supported by actual run times. A reference implementation is provided together with instructions on how to reproduce the main results of this paper. Reproducibility of computational results. This paper has been awarded the “SIAM Reproducibility Badge: Code and data available” as recognition that the authors have followed reproducibility principles valued by SISC and the scientific computing community. Code and data that allow readers to reproduce the results in this paper are available at https://github.com/un110076/DifferentialInversion and in the supplementary materials (DifferentialInversion-main.zip [12.5KB]).

Estimate of Koopman Modes and Eigenvalues with Kalman Filter

Ningxin Liu — 2026-03-10T07:00:00Z

SIAM Journal on Scientific Computing, Volume 48, Issue 2, Page A512-A539, April 2026.
Abstract. Dynamic mode decomposition (DMD) is a data-driven method of extracting spatial-temporal coherent modes from complex systems and providing an equation-free architecture to model and predict systems. However, in practical applications, the accuracy of DMD can be limited in extracting dynamical features due to sensor noise in measurements. We develop an adaptive method to constantly update dynamic modes and eigenvalues from noisy measurements arising from discrete systems. Our method is based on the Ensemble Kalman filter owing to its capability of handling time-varying systems and nonlinear observables. Our method can be extended to nonautonomous dynamical systems, accurately recovering short-time eigenvalue-eigenvector pairs and observables. Theoretical analysis shows that the estimation is accurate in long term data misfit. We demonstrate the method on both autonomous and nonautonomous dynamical systems to show its effectiveness.

High-Order Finite Element Methods for Three-Dimensional Multicomponent Convection-Diffusion

Aaron Baier-Reinio — 2026-03-13T07:00:00Z

SIAM Journal on Scientific Computing, Volume 48, Issue 2, Page A540-A567, April 2026.
Abstract. We derive and analyze a broad class of finite element methods for numerically simulating the stationary, low Reynolds number flow of concentrated mixtures of several distinct chemical species in a common thermodynamic phase. The underlying partial differential equations that we discretize are the Stokes–Onsager–Stefan–Maxwell (SOSM) equations, which model bulk momentum transport and multicomponent diffusion within ideal and nonideal mixtures. Unlike previous approaches, the methods are straightforward to implement in two and three spatial dimensions and allow for high-order finite element spaces to be employed. The key idea in deriving the discretization is to suitably reformulate the SOSM equations in terms of the species mass fluxes and chemical potentials and discretize these unknown fields using stable [math]–[math] finite element pairs. We prove that the methods are convergent and yield a symmetric linear system for a Picard linearization of the SOSM equations, which staggers the updates for concentrations and chemical potentials. We also discuss how the proposed approach can be extended to the Newton linearization of the SOSM equations, which requires the simultaneous solution of mole fractions, chemical potentials, and other variables. Our theoretical results are supported by numerical experiments, and we present an example of a physical application involving the microfluidic nonideal mixing of hydrocarbons. Reproducibility of computational results. This paper has been awarded the “SIAM Reproducibility Badge: Code and data available” as a recognition that the authors have followed reproducibility principles valued by SISC and the scientific computing community. Code and data that allow readers to reproduce the results in this paper are available at https://zenodo.org/records/16416181.

Interior Penalty Discontinuous Galerkin Methods for the Nearly Incompressible Elasticity Eigenvalue Problem with Heterogeneous Media

Arbaz Khan — 2026-03-19T07:00:00Z

SIAM Journal on Scientific Computing, Volume 48, Issue 2, Page A568-A599, April 2026.
Abstract. This paper studies the family of interior penalty discontinuous Galerkin methods for solving the Herrmann formulation of the linear elasticity eigenvalue problem in heterogeneous media. By employing a weighted Lamé coefficient norm within the framework of noncompact operators theory, we prove convergence of both continuous and discrete eigenvalue problems as the mesh size approaches zero, independently of the Lamé constants. Additionally, we conduct an a posteriori analysis and propose a reliable and efficient estimator. Our theoretical findings are supported by numerical experiments.

A Preconditioned Riemannian Conjugate Gradient Method to Compute Ground States of Spin-1 Bose–Einstein Condensates

Wen Huang — 2026-04-01T07:00:00Z

SIAM Journal on Scientific Computing, Volume 48, Issue 2, Page A600-A623, April 2026.
Abstract. In this paper, we propose a preconditioned Riemannian conjugate gradient method for computing ground states of spin-1 Bose-Einstein condensates, which can be reformulated as an optimization problem with the mass and magnetization constraints. The energy functional and constraints can be discretized using Fourier pseudospectral schemes, thereby transforming the problem into an optimization task on a manifold. We derive vector transports by differentiating three existing retractions and propose an initial step size selection strategy based on the second-order approximation of the energy function for the proposed Riemannian optimization algorithm. In addition, a preconditioner is derived and used. To further accelerate the convergence, we combine the proposed algorithm with a multigrid method. Numerical experiments demonstrate the efficiency and accuracy of the proposed method while confirming the usefulness of our preconditioning technique and the effectiveness of the proposed strategy for initial step size selection.

Fast Measure Modification of Orthogonal Polynomials via Matrices with Displacement or Hierarchical Off-Diagonal Low-Rank Structure

Karim Gumerov — 2026-04-02T07:00:00Z

SIAM Journal on Scientific Computing, Volume 48, Issue 2, Page A624-A645, April 2026.
Abstract. It is well known that matrices with low Hessenberg-structured displacement rank enjoy fast algorithms for certain matrix factorizations. We show how [math] principal finite sections of the Gram matrix for the orthogonal polynomial measure modification problem have such a displacement structure, unlocking a collection of fast algorithms for computing connection coefficients (as the upper-triangular Cholesky factor) between a known orthogonal polynomial family and the modified family. In general, the [math] complexity is reduced to [math], and if the symmetric Gram matrix has upper and lower bandwidth [math], then the [math] complexity for a banded Cholesky factorization is reduced to [math]. In the case of modified Chebyshev polynomials, we show that the Gram matrix is a symmetric Toeplitz-plus-Hankel matrix, and if the modified Chebyshev moments decay algebraically, then a hierarchical off-diagonal low-rank structure is observed in the Gram matrix, enabling a further reduction in the complexity of an approximate Cholesky factorization powered by randomized numerical linear algebra.

A BDDC Preconditioner for the Cardiac EMI Model in Three Dimensions

Fritz Göbel — 2026-04-03T07:00:00Z

SIAM Journal on Scientific Computing, Volume 48, Issue 2, Page A646-A667, April 2026.
Abstract. We analyze a balancing domain decomposition by constraints (BDDC) preconditioner for the solution of three-dimensional composite discontinuous Galerkin discretizations of reaction-diffusion systems of ordinary and partial differential equations arising in cardiac cell-by-cell models like the extracellular space, membrane and intracellular space (EMI) model. These microscopic models are essential for the understanding of events in aging and structurally diseased hearts, which macroscopic models relying on homogenized descriptions of the cardiac tissue, like monodomain and bidomain models, fail to adequately represent. The modeling of each individual cardiac cell results in discontinuous global solutions across cell boundaries, requiring the careful construction of dual and primal spaces for the BDDC preconditioner. We provide a scalable condition number bound for the precondition operator and validate the theoretical results with extensive numerical experiments.

Radial Basis Function Techniques for Neural Field Models on Surfaces

Sage B. Shaw — 2026-04-03T07:00:00Z

SIAM Journal on Scientific Computing, Volume 48, Issue 2, Page A668-A694, April 2026.
Abstract. We present a numerical framework for solving neural field equations on surfaces using radial basis function (RBF) interpolation and quadrature. Neural field models describe the evolution of macroscopic brain activity, but modeling studies often overlook the complex geometry of curved cortical domains. Traditional numerical methods, such as finite element or spectral methods, can be computationally expensive and challenging to implement on irregular domains. In contrast, RBF-based methods provide a flexible alternative by offering interpolation and quadrature schemes that efficiently handle arbitrary geometries with high-order accuracy. We first develop an RBF-based interpolatory projection framework for neural field models on general smooth surfaces. Quadratures for both flat and curved domains are derived in detail, ensuring high-order accuracy and stability as they depend on RBF hyper-parameters (basis functions, augmenting polynomials, and stencil size). Through numerical experiments, we demonstrate the convergence of our method, highlighting its advantages over traditional approaches in terms of flexibility and accuracy. We conclude with an exposition of numerical simulations of spatiotemporal activity on complex surfaces, illustrating the method’s ability to capture complex wave propagation patterns. Reproducibility of computational results. This paper has been awarded the “SIAM Reproducibility Badge: Code and data available” as a recognition that the authors have followed reproducibility principles valued by SISC and the scientific computing community. Code and data that allow readers to reproduce the results in this paper are available at https://github.com/shawsa/neural-fields-rbf and in the supplementary materials (M175378_Supplementary_Materials_1.pdf [657KB]).

Oscillation-Eliminating Central DG Schemes on Cartesian Meshes for Hyperbolic Conservation Laws

Manting Peng — 2026-04-06T07:00:00Z

SIAM Journal on Scientific Computing, Volume 48, Issue 2, Page A695-A725, April 2026.
Abstract. This paper proposes and analyzes the OECDG method, a novel central discontinuous Galerkin (CDG) scheme that integrates the strengths of the oscillation-eliminating (OE) approach [M. Peng, Z. Sun, and K. Wu, Math. Comp., 94 (2025), pp. 1147–1198] within the CDG framework for general hyperbolic conservation laws. The OECDG method incorporates a new OE procedure, along with an innovative dual damping mechanism, to enhance oscillation control in CDG schemes. Unlike OE procedures in standard DG methods that use intercell jumps as a modal filter, our new OE mechanism draws inspiration from CDG-based numerical dissipation and employs a convex blending strategy, leveraging overlapping solutions to enhance stencil compactness and resolution without characteristic decomposition. We rigorously derive optimal error estimates for the fully discrete OECDG method through several key theoretical advancements, filling a gap in the error analysis of fully discrete CDG schemes, including original linear CDG methods without oscillation control. First, we establish the approximate skew-symmetry and weak boundedness of the CDG spatial discretization, which underpins the fully discrete stability analysis. Utilizing these properties, we then prove the linear stability of CDG methods coupled with Runge–Kutta time discretization via matrix transfer techniques. These foundational results enable us to derive fully discrete optimal error estimates for the OECDG method—a challenging task due to the method’s nonlinear nature, even for linear advection equations. Extensive numerical experiments validate the theoretical findings and demonstrate the efficacy of the OECDG method across a variety of hyperbolic conservation laws, including linear and nonlinear problems such as the convection equation, the Burgers equation, a traffic flow model, and Euler equations. The results confirm the OECDG method’s ability to achieve optimal convergence rates, robustly eliminate spurious oscillations, and accurately capture complex wave structures across diverse test cases.

A Scaling and Recovering Algorithm for the Matrix [math]-Functions

Awad H. Al-Mohy — 2026-04-06T07:00:00Z

SIAM Journal on Scientific Computing, Volume 48, Issue 2, Page A726-A747, April 2026.
Abstract. A new scaling and recovering algorithm is proposed for simultaneously computing the matrix [math]-functions that arise in exponential integrator methods for the numerical solution of certain first-order systems of ordinary differential equations. The algorithm initially scales the input matrix down by a nonnegative integer power of two, and then evaluates the [math] diagonal Padé approximant to [math], where [math] is the largest index of interest. The remaining [math] Padé approximants to [math], [math], are obtained implicitly via a recurrence relation. The effect of scaling is subsequently recovered using the double-argument formula. A rigorous backward error analysis, based on the [math] Padé approximant to the exponential, enables sharp bounds on the relative backward errors. These bounds are expressed in terms of the sequence [math], which can be much smaller than [math] for nonnormal matrices. The scaling parameter and the degrees of the Padé approximants are selected to minimize the overall computational cost, which benefits from the sharp bounds and the optimal evaluation schemes for diagonal Padé approximants. Furthermore, if the input matrix is (quasi-)triangular, the algorithm exploits its structure in the recovering phase. Numerical experiments demonstrate the superiority of the proposed algorithm over existing alternatives in both accuracy and efficiency. Reproducibility of computational results. This paper has been awarded the “SIAM Reproducibility Badge: Code and data available” as a recognition that the authors have followed reproducibility principles valued by SISC and the scientific computing community. Code and data that allow readers to reproduce the results in this paper are available at https://github.com/xiaobo-liu/phi_funm and in the supplementary materials (phi_funm-main.zip [8.39MB]).

A Simple and Fast Finite Difference Method for the Integral Fractional Laplacian of Variable Order

Zhaopeng Hao — 2026-04-06T07:00:00Z

SIAM Journal on Scientific Computing, Volume 48, Issue 2, Page A748-A783, April 2026.
Abstract. The numerical evaluation of multidimensional variable-order fractional Laplacians presents significant challenges due to the inherent singularity of the defining integral. We develop a simple yet effective finite difference scheme for approximating this hypersingular integral operator. Our approach achieves second-order convergence while maintaining straightforward implementation. For the associated fractional Poisson equation with variable-order Laplacian, we first establish solution wellposedness before constructing a corresponding finite difference discretization based on our second-order approximation. A fast algorithm implementation is presented to enhance computational efficiency. Numerical experiments across multiple dimensions validate our theoretical results, demonstrating both the accuracy of our second-order convergence and the practical efficiency of the proposed method.

A Nonuniform Fast Hankel Transform

Paul G. Beckman — 2026-04-06T07:00:00Z

SIAM Journal on Scientific Computing, Volume 48, Issue 2, Page A784-A803, April 2026.
Abstract. We describe a fast algorithm for computing discrete Hankel transforms of moderate orders from [math] nonuniform points to [math] nonuniform frequencies in [math] operations. Our approach combines local and asymptotic Bessel function expansions with nonuniform fast Fourier transforms. The order of each expansion is adjusted automatically according to error analysis to obtain any desired precision [math]. Several numerical examples are provided which demonstrate the speed and accuracy of the algorithm in multiple regimes and applications. Reproducibility of computational results. This paper has been awarded the “SIAM Reproducibility Badge: Code and data available” as a recognition that the authors have followed reproducibility principles valued by SISC and the scientific computing community. Code and data that allow readers to reproduce the results in this paper are available at https://github.com/pbeckman/FastHankelTransform.jl and in the supplementary materials (FastHankelTransform-jl-main.zip [21.4KB]).

A Time-Dependent Phase Space Filter for Anisotropic Wave Equations on Unbounded Domains

Thomas Hagstrom — 2026-04-07T07:00:00Z

SIAM Journal on Scientific Computing, Volume 48, Issue 2, Page A804-A827, April 2026.
Abstract. For wave propagation problems in homogeneous, isotropic media, exponentially convergent domain truncation techniques based on perfectly matched damping layers or optimal rational approximations to the exact radiation operator are known. However, these methods fail for systems with wave families satisfying general dispersion relations, forcing practitioners to resort to ad hoc procedures based on grid stretching and artificial damping. Here we propose a new method for constructing convergent approximations on truncated domains, the phase space filter, which unlike other methods is completely general and mathematically justified. Based on the fact that outgoing waves can be characterized as waves located near the boundary of the computational domain with group velocities pointing outward, the key idea of the phase space filtering algorithm consists of applying a filter to the solution that removes outgoing waves only. The method introduced in this work is a simplified version of the original phase space filter, which was proposed for the Schrödinger equation. The method is applied to anisotropic wave models for which existing techniques are unstable, namely free-space problems governed in the far field by the Euler equations linearized about a uniform mean flow and Maxwell’s equations in an anisotropic medium. We also present experiments in waveguide geometry and for isotropic problems which indicate that a more sophisticated multiscale extension of the algorithm is needed to obtain good long-time accuracy for waveguide problems, or to match the performance of the optimized local radiation conditions available in the isotropic case. Theoretical results concerning the convergence and computational costs of the phase space filter are discussed and stability is proven. Reproducibility of computational results. This paper has been awarded the “SIAM Reproducibility Badge: Code and data available” as a recognition that the authors have followed reproducibility principles valued by SISC and the scientific computing community. Code and data that allow readers to reproduce the results in this paper are available at https://zenodo.org/records/17237326.

Algebraic Multigrid with Filtering: An Efficient Preconditioner for Interior Point Methods in Large-Scale Contact Mechanics Optimization

Socratis Petrides — 2026-04-08T07:00:00Z

SIAM Journal on Scientific Computing, Volume 48, Issue 2, Page A828-A852, April 2026.
Abstract. Large-scale contact mechanics simulations are crucial in many engineering fields such as structural design and manufacturing. In the frictionless case, contact can be modeled by minimizing an energy functional; however, these problems are often nonlinear, nonconvex, and increasingly difficult to solve as mesh resolution increases. In this work, we employ a Newton-based interior-point (IP) filter line-search method, an effective approach for large-scale constrained optimization. While this method converges rapidly, each iteration requires solving a large saddle-point linear system that becomes ill-conditioned as the optimization process converges, largely due to IP treatment of the contact constraints. Such ill-conditioning can hinder solver scalability and increase iteration counts with mesh refinement. To address this, we introduce a novel preconditioner, algebraic multigrid with filtering (AMGF), tailored to the Schur complement of the saddle-point system. Building on the classical AMG solver, commonly used for elasticity, we augment it with a specialized subspace correction that filters near null space components introduced by contact interface constraints. Through theoretical analysis and numerical experiments on a range of linear and nonlinear contact problems, we demonstrate that the proposed solver achieves mesh independent convergence and maintains robustness against the ill-conditioning that notoriously plagues IP methods. These results indicate that AMGF makes contact mechanics simulations more tractable and broadens the applicability of Newton-based IP methods in challenging engineering scenarios. More broadly, AMGF is well suited for problems, optimization or otherwise, where solver performance is limited by a low-dimensional subspace, such as those arising from localized constraints, interface conditions, or model heterogeneities. This makes the method widely applicable beyond contact mechanics and constrained optimization.

Localized Dynamic Mode Decomposition with Temporally Adaptive Segmentation

Qiuqi Li — 2026-04-08T07:00:00Z

SIAM Journal on Scientific Computing, Volume 48, Issue 2, Page A853-A877, April 2026.
Abstract. Dynamic mode decomposition (DMD) is a widely used data-driven method for predicting the future states of dynamical systems. However, its standard formulation often struggles with poor long-term predictive accuracy. To address this limitation, we propose a localized DMD (LDMD) framework that improves prediction performance by integrating DMD’s strong linear forecasting capabilities with time-domain segmentation techniques. In this framework, the temporal domain is segmented into multiple sub-intervals, within which snapshot matrices are constructed and localized predictions are performed. We first present the localized DMD method with predefined segmentation, and then explore an adaptive segmentation strategy to further enhance computational efficiency and prediction robustness. Furthermore, we conduct an error analysis that provides the upper bound of the local and global truncation error for the proposed framework. The effectiveness of LDMD is demonstrated on four benchmark problems—Burgers’, Allen–Cahn, nonlinear Schrödinger, and Maxwell’s equations. Numerical results show that LDMD significantly enhances long-term predictive accuracy while preserving high computational efficiency.

Parameter-Robust Preconditioner for Stokes-Darcy Coupled Problem with Lagrange Multiplier

Xiaozhe Hu — 2026-04-09T07:00:00Z

SIAM Journal on Scientific Computing, Volume 48, Issue 2, Page A878-A904, April 2026.
Abstract. In this paper, we propose a parameter-robust preconditioner for the coupled Stokes–Darcy problem equipped with various boundary conditions, enforcing the mass conservation at the interface via a Lagrange multiplier. We rigorously establish that the coupled system is well-posed with respect to physical parameters and mesh size and provides a framework for constructing parameter-robust preconditioners. Furthermore, we analyze the convergence behavior of the minimal residual method in the presence of small outlier eigenvalues linked to specific boundary conditions, which can lead to slow convergence or stagnation. To address this issue, we employ deflation techniques to accelerate the convergence. Finally, numerical experiments confirm the effectiveness and robustness of the proposed approach.

Bayesian D-Optimal Experimental Designs via Column Subset Selection

Srinivas Eswar — 2026-04-09T07:00:00Z

SIAM Journal on Scientific Computing, Volume 48, Issue 2, Page A905-A928, April 2026.
Abstract. This paper tackles optimal sensor placement for Bayesian linear inverse problems, a popular version of the more general optimal experimental design (OED) problem, using the D-optimality criterion. This is done by establishing connections between sensor placement and the column subset selection problem (CSSP), which is a well-studied problem in numerical linear algebra (NLA). In particular, we use the Golub–Klema–Stewart (GKS) approach, which involves computing the truncated singular value decomposition (SVD) followed by a pivoted QR factorization on the right singular vectors. The algorithms are further accelerated by using randomization to compute the low-rank approximation as well as for sampling the indices. The resulting algorithms are robust, computationally efficient, amenable to parallelization, require virtually no parameter tuning, and come with strong theoretical guarantees. One of the proposed algorithms is also adjoint-free, which is beneficial in situations where the adjoint is expensive to evaluate or is not available. Additionally, we develop a method for data completion without solving the inverse problem. Numerical experiments on model inverse problems involving the heat equation and seismic tomography in two spatial dimensions demonstrate the performance of our approaches. Reproducibility of computational results. This paper has been awarded the “SIAM Reproducibility Badge: Code and Data Available” as a recognition that the authors have followed reproducibility principles valued by SISC and the scientific computing community. Code and data that allow readers to reproduce the results in this paper are available at https://github.com/RandomizedOED/css4oed and in the supplementary materials (M164986_Supplementary_Materials_2.pdf [455KB], css4oed-main.zip [820KB]).

Truncated Huber Penalty for Sparse Signal Recovery with Convergence Analysis

Li Yang — 2026-04-10T07:00:00Z

SIAM Journal on Scientific Computing, Volume 48, Issue 2, Page A929-A957, April 2026.
Abstract. Sparse signal recovery from underdetermined systems presents significant challenges when using conventional [math] and [math] penalties, primarily due to computational complexity and estimation bias. This paper introduces a truncated Huber penalty, a nonconvex metric that effectively bridges the gap between unbiased sparse recovery and differentiable optimization. The proposed penalty applies quadratic regularization to small entries while truncating large magnitudes, avoiding nondifferentiable points at optimal solutions. Theoretical analysis demonstrates that, for an appropriately chosen threshold, any [math]-sparse solution recoverable via conventional penalties remains a local optimum under the truncated Huber function. This property allows the exact and robust recovery theories developed for other penalty regularization functions to be directly extended to the truncated Huber function. To solve the optimization problem, we develop a block coordinate descent (BCD) algorithm with finite-step convergence guarantees under spark conditions. Numerical experiments validate the effectiveness and robustness of the proposed method in both synthetic and real scenarios. Furthermore, we demonstrate the flexibility of the truncated Huber framework through two extensions: one to an adaptively weighted variant inspired by sorted penalties, and the other to the gradient domain for applications such as signal denoising and image smoothing.

Integral Kernel Methods for Nonlinear Parabolic-Elliptic Systems

Henrique B. N. Monteiro — 2026-04-13T07:00:00Z

SIAM Journal on Scientific Computing, Volume 48, Issue 2, Page A958-A983, April 2026.
Abstract. Nonlinear parabolic-elliptic systems arise in many physical, biological, and chemical phenomena such as chemotaxis, ion transport, self-gravitating particles, and Brownian vortices. Existing methods struggle with the strong coupling and high nonlinearity and nonlocality of some of these systems, especially the ill-conditioned, convection-dominated problems. To overcome numerical difficulties, current approaches rely on initial guesses, preconditioning, or iterative techniques with no convergence guarantees. They might suffer from poor scalability, large memory usage, and difficulty to parallelize. Inspired by the connection of parabolic-elliptic systems to stochastic processes, we introduce a novel meshless, monolithic, and fully explicit method that naturally encapsulates the elliptic and parabolic operators into a single step which updates each node deterministically with global information. By being fully quadrature-based, it avoids solving systems of discretized equations and does not utilize initial guesses or preconditioning, while requiring little memory and being easy to parallelize. We first derive the method in an integral kernel formulation with quadratic complexity in the number of integration nodes and then leverage kernel-independent fast multipole methods (FMM) to present a scalable algorithm with linear complexity. We provide numerical examples for the Poisson–Nernst–Planck equations in one, two, and three dimensions, together with the derivation of the integral kernel for each case. The examples demonstrate the fast convergence and scalability of the FMM-accelerated algorithm, as well as its suitability for convection-dominated problems, making it competitive against traditional PDE solvers.

A Space-Time Adaptive Boundary Element Method for the Wave Equation

A. Aimi — 2026-04-13T07:00:00Z

SIAM Journal on Scientific Computing, Volume 48, Issue 2, Page A984-A1005, April 2026.
Abstract. This article initiates the study of space-time adaptive mesh refinements for time-dependent boundary element formulations of wave equations. Based on error indicators of residual type, we formulate an adaptive boundary element procedure for acoustic soft-scattering problems with local tensor-product refinements of the space-time mesh. We discuss the algorithmic challenges and investigate the proposed method in numerical experiments. In particular, we study the performance and improved convergence rates with respect to the energy norm for problems dominated by spatial, temporal, or traveling singularities of the solution. The efficiency of the considered rigorous and heuristic a posteriori error indicators is discussed.

A Hierarchical Approach for Multicontinuum Homogenization in High Contrast Media

Wei Xie — 2026-04-16T07:00:00Z

SIAM Journal on Scientific Computing, Volume 48, Issue 2, Page A1006-A1027, April 2026.
Abstract. A recently developed upscaling technique, the multicontinuum homogenization method, has gained significant attention for its effectiveness in modeling complex multiscale systems. This method defines multiple continua based on distinct physical properties and solves a series of constrained cell problems to capture localized information for each continuum. However, solving all these cell problems on very fine grids at every macroscopic point is computationally expensive, which is a common limitation of most homogenization approaches for nonperiodic problems. To address this challenge, we propose a hierarchical multicontinuum homogenization framework. The core idea is to define hierarchical macroscopic points and solve the constrained problems on grids of varying resolutions. The local solutions are decomposed into the linear interpolation of contributions inherited from preceding levels and an additional correction term. This combination is substituted into the original constrained problems, and the correction term is resolved using finite element (FE) grids of varying sizes depending on the level of the macropoint. By normalizing the computational cost of fully resolving the local problem to [math], we establish that our approach incurs a cost of [math], highlighting substantial computational savings across hierarchical layers [math], coarsening factor [math], and spatial dimension [math]. Numerical experiments validate the effectiveness of the proposed method in media with slowly varying properties, underscoring its potential for efficient multiscale modeling. Reproducibility of computational results. This paper has been awarded the “SIAM Reproducibility Badge: Code and data available” as recognition that the authors have followed reproducibility principles valued by SISC and the scientific computing community. Code and data that allow readers to reproduce the results in this paper are available at https://github.com/xieweidc/hierarchical_mh_elliptic.git and in the supplementary materials (hierarchical_mh_elliptic-main.zip [19.8KB]).

A Finite Element Method for Maxwell’s Transmission Eigenvalue Problem in Anisotropic Media

Jiayu Han — 2026-04-17T07:00:00Z

SIAM Journal on Scientific Computing, Volume 48, Issue 2, Page A1028-A1049, April 2026.
Abstract. In this paper, we introduce a finite element method employing the Nedéléc element space for solving Maxwell’s transmission eigenvalue problem in anisotropic media. The well-posedness of the source problems is derived using the [math]-coercivity approach. We discuss the discrete compactness property of the finite element space under the case of anisotropic coefficients and conduct a finite element error analysis for the proposed approach. Additionally, we present some numerical examples to support the theoretical result. Reproducibility of computational results. This paper has been awarded the “SIAM Reproducibility Badge: Code and data available” as recognition that the authors have followed reproducibility principles valued by SISC and the scientific computing community. Code and data that allow readers to reproduce the results in this paper are available at https://github.com/hanjiayu126/MTEP_bd and in the supplementary material (MTEP_bd-main.zip [15.2KB]).

A Fast and Accurate Solver for the Fractional Fokker–Planck Equation with Dirac-Delta Initial Conditions

Qihao Ye — 2026-04-20T07:00:00Z

SIAM Journal on Scientific Computing, Volume 48, Issue 2, Page A1050-A1074, April 2026.
Abstract. The classical Fokker–Planck equation (FPE) is a key tool in physics for describing systems influenced by drag forces and Gaussian noise, with applications spanning multiple fields. We consider the fractional Fokker–Planck equation (FFPE), which models the time evolution of probability densities for systems driven by Lévy processes, relevant in scenarios where Gaussian assumptions fail. The paper presents an efficient and accurate numerical approach for the free-space FFPE with constant coefficients and Dirac-delta initial conditions. This method utilizes the integral representation of the solutions and enables the efficient handling of very high-dimensional problems using fast algorithms. Our work is the first to present a high-precision numerical solver for the free-space FFPE with Dirac-delta initial conditions. This opens the door for future research on more complex scenarios, including those with variable coefficients and other types of initial conditions. Reproducibility of computational results. This paper has been awarded the “SIAM Reproducibility Badge: Code and data available” as a recognition that the authors have followed reproducibility principles valued by SISC and the scientific computing community. Code and data that allow readers to reproduce the results in this paper are available at https://github.com/ACMathX/FFPEDDIC and in the supplementary materials (M168290_Supplementary_Materials.pdf [675KB], FFPEDDIC-main.zip [3.68MB]).

Shape Optimization with Ventcel Transmission Condition: Application to the Design of a Heat Exchanger

Fabien Caubet — 2026-03-04T08:00:00Z

SIAM Journal on Scientific Computing, Volume 48, Issue 2, Page B113-B140, April 2026.
Abstract. This paper aims to optimize the shape of a fluid-to-fluid heat exchanger in order to maximize heat exchange under constraints of energy dissipation and volume. The novelty consists in taking into account the thin layer separating the two fluids by using Ventcel-type second-order transmission conditions. The physical model is then a weakly coupled problem between the steady-state Navier–Stokes equations for the dynamics of the two fluids dynamics and the convection-diffusion equation for the heat. We provide a shape sensitivity analysis and characterize the shape derivatives involved. Finally, we demonstrate the feasibility and effectiveness of the proposed method through three-dimensional numerical simulations.

Exploiting Inexact Computations in Multilevel Monte Carlo and Other Sampling Methods

Josef Martínek — 2026-03-05T08:00:00Z

SIAM Journal on Scientific Computing, Volume 48, Issue 2, Page B141-B164, April 2026.
Abstract. Multilevel sampling methods, such as multilevel and multifidelity Monte Carlo, multilevel stochastic collocation, or delayed acceptance Markov chain Monte Carlo, have become standard uncertainty quantification (UQ) tools for a wide class of forward and inverse problems. The underlying idea is to achieve faster convergence by leveraging a hierarchy of models, such as partial differential equation (PDE) or stochastic differential equation (SDE) discretisations with increasing accuracy. By optimally redistributing work among the levels, multilevel methods can achieve significant performance improvement compared to single level methods working with one high-fidelity model. Intuitively, approximate solutions on coarser levels can tolerate large computational error without affecting the overall accuracy. We show how this can be used in high-performance computing applications to obtain a significant performance gain. As a use case, we analyze the computational error in the standard multilevel Monte Carlo method and formulate an adaptive algorithm which determines a minimum required computational accuracy on each level of discretisation. We show two examples of how the inexactness can be converted into actual gains using an elliptic PDE with lognormal random coefficients. Using a low precision sparse direct solver combined with iterative refinement results in a simulated gain in memory references of up to [math] compared to the reference double precision solver; while using a MINRES iterative solver, a practical speedup of up to [math] in terms of FLOPs is achieved. These results provide a step in the direction of energy-aware scientific computing, with significant potential for energy savings. Reproducibility of computational results. This paper has been awarded the “SIAM Reproducibility Badge: Code and data available” as a recognition that the authors have followed reproducibility principles valued by SISC and the scientific computing community. Code and data that allow readers to reproduce the results in this paper are available at https://github.com/josef-martinek/mpml and in the supplementary materials (mpml-main.zip [72.3KB]).

Helicity-Preserving Finite Element Discretization for Magnetic Relaxation

Mingdong He — 2026-03-10T07:00:00Z

SIAM Journal on Scientific Computing, Volume 48, Issue 2, Page B165-B183, April 2026.
Abstract. The Parker conjecture, which explores whether magnetic fields in perfectly conducting plasmas can develop tangential discontinuities during magnetic relaxation, remains an open question in astrophysics. Helicity conservation provides a topological barrier during relaxation, preventing topologically nontrivial initial data relaxing to trivial solutions; preserving this mechanism discretely over long time periods is therefore crucial for numerical simulation. This work presents an energy- and helicity-preserving finite element discretization for the magneto-frictional system for investigating the Parker conjecture. The algorithm preserves a discrete version of the topological barrier and a discrete Arnold inequality. We also propose extensions of the notion of helicity and the Arnold inequality to certain kinds of topologically nontrivial domains. Numerical experiments demonstrate that helicity preservation is crucial in obtaining physically meaningful simulations of magnetic relaxation, providing an example where structure-preserving schemes are necessary. Reproducibility of computational results. This paper has been awarded the “SIAM Reproducibility Badge: Code and data available” as a recognition that the authors have followed reproducibility principles valued by SISC and the scientific computing community. Code and data that allow readers to reproduce the results in this paper are available at https://zenodo.org/records/15302724 and https://zenodo.org/records/16797562.

Coupled Eikonal Problems to Model Cardiac Reentries in Purkinje Network and Myocardium

Samuele Brunati — 2026-03-16T07:00:00Z

SIAM Journal on Scientific Computing, Volume 48, Issue 2, Page B184-B207, April 2026.
Abstract. We propose a novel partitioned scheme based on Eikonal equations to model the coupled propagation of the electrical signal in the His–Purkinje system and in the myocardium for cardiac electrophysiology. This scheme allows, for the first time in Eikonal-based modeling, to capture all possible signal reentries between the Purkinje network and the cardiac muscle that may occur under pathological conditions. As part of the proposed scheme, we introduce a new pseudo-time method for the Eikonal-diffusion problem in the myocardium, to correctly enforce electrical stimuli coming from the Purkinje network. We test our approach by performing numerical simulations of cardiac electrophysiology in a real biventricular geometry, under both pathological and therapeutic conditions, to demonstrate its flexibility, robustness, and predictive capabilities.

Fast and Accurate Intersections on a Sphere

Hongyu Chen — 2026-04-10T07:00:00Z

SIAM Journal on Scientific Computing, Volume 48, Issue 2, Page B208-B232, April 2026.
Abstract. We introduce a fast, high-precision algorithm for calculating intersections between great circle arcs and lines of constant latitude on the unit sphere. We first propose a simplified intersection point formula with improved speed and numerical robustness over the ones traditionally implemented in geoscience software. We then show how algorithms based on the concept of error-free transformations (EFT) can be applied to evaluate this formula within a relative error bound that is on the order of machine precision. We demonstrate that, with a vectorized and parallelized implementation, this enhanced accuracy is achieved with no compute time overhead compared to a direct calculation in hardware floating point, making our algorithm suitable for performance-sensitive applications like regridding of high-resolution climate data. In contrast, evaluating our formula using high-precision data types like quadruple precision and arbitrary precision, or using the robust intersection computation routines from the Computational Geometry Algorithms Library, leads to significant computational overhead, especially since these alternatives inhibit vectorization. More generally, our work demonstrates how EFT techniques can be combined and extended to implement nontrivial geometric calculations with high accuracy and speed. Reproducibility of computational results. This paper has been awarded the “SIAM Reproducibility Badge: Code and data available” as a recognition that the authors have followed reproducibility principles valued by SISC and the scientific computing community. Code and data that allow readers to reproduce the results in this paper are available at https://github.com/hongyuchen1030/accurate-operator-reproducibility/tree/main and in the supplementary materials (M173761_Supplementary_Materials_1.pdf [265KB], accurate-operator-reproducibility-main.zip [199KB]).

Efficient, Decoupled, Second and Third-Order Implicit-Explicit Runge–Kutta Schemes with Original-Form Energy Stability for the Anisotropic Phase-Field Dendritic Growth Model

Jianan Li — 2026-04-15T07:00:00Z

SIAM Journal on Scientific Computing, Volume 48, Issue 2, Page B233-B261, April 2026.
Abstract. We consider numerical approximations of the anisotropic phase-field dendritic crystal growth model, a highly nonlinear system comprising the anisotropic Allen–Cahn equation and the thermal equation, coupled through the latent heat effect. Despite various efforts to develop efficient numerical approaches, significant challenges remain in constructing high-order (beyond second-order) and fully decoupled numerical schemes with provable energy stability, particularly with respect to the energy form corresponding to the original system. To address these challenges, we have developed a novel, linear, high-order (second or third) fully decoupled implicit-explicit Runge–Kutta scheme. Our schemes preserve energy stability in their original form, with only minor Lipschitiz modifications to the nonlinear potentials, which differ fundamentally from approaches that reformulate the entire free energy via auxiliary variables. The key to achieving energy stability lies in utilizing stabilization technique and various coefficient matrices for different nonlinear terms. We propose specific requirements for these coefficient matrices and introduce a search strategy to identify appropriate schemes, rigorously proving their energy stability. Numerous 2D and 3D simulations have been performed to demonstrate the stability and accuracy of the proposed schemes, numerically.

Robust Optimal Experimental Design of Infinite-Dimensional Bayesian Nonlinear Inverse Problems

Abhijit Chowdhary — 2026-04-20T07:00:00Z

SIAM Journal on Scientific Computing, Volume 48, Issue 2, Page B262-B288, April 2026.
Abstract. We consider robust optimal experimental design (ROED) for nonlinear Bayesian inverse problems governed by partial differential equations (PDEs). An optimal design is one that maximizes some utility quantifying the quality of the solution of an inverse problem. However, the optimal design is dependent on elements of the inverse problem such as the simulation model, the prior, or the measurement error model. ROED aims to produce an optimal design that is aware of the additional uncertainties encoded in the inverse problem and remains optimal even after variations in them. We follow a worst-case scenario approach to develop a new framework for robust optimal design of nonlinear Bayesian inverse problems. The proposed framework (a) is scalable and designed for infinite-dimensional Bayesian nonlinear inverse problems constrained by PDEs; (b) develops efficient approximations of the utility, namely the expected information gain; (c) employs eigenvalue sensitivity techniques to develop analytical forms and efficient evaluation methods of the gradient of the utility with respect to the uncertainties against which we wish to be robust; and (d) employs a probabilistic optimization paradigm that properly defines and efficiently solves the resulting combinatorial max-min optimization problem. The effectiveness of the proposed approach is illustrated for optimal sensor placement problem in an inverse problem governed by an elliptic PDE.

Discovery of Probabilistic Dirichlet-to-Neumann Maps on Graphs

Adrienne M. Propp — 2026-03-20T07:00:00Z

SIAM Journal on Scientific Computing, Volume 48, Issue 2, Page C191-C215, April 2026.
Abstract. Dirichlet-to-Neumann maps enable the coupling of multiphysics simulations across computational subdomains by ensuring continuity of state variables and fluxes at artificial interfaces. We present a novel method for learning Dirichlet-to-Neumann maps on graphs using Gaussian processes, specifically for problems where the data obey a conservation law arising from an underlying partial differential equation. Our approach combines discrete exterior calculus and nonlinear optimal recovery to infer relationships between vertex and edge values. This framework yields data-driven predictions with uncertainty quantification across the entire graph, even when observations are limited to a subset of vertices and edges. By minimizing the reproducing kernel Hilbert space norm while penalizing kernel complexity through maximum likelihood estimation, our method ensures that the resulting surrogate strictly enforces conservation laws without overfitting. We demonstrate our method on two representative applications: subsurface flow in fracture networks and arterial blood flow. The results demonstrate that the method maintains high accuracy and well-calibrated uncertainty estimates even under severe data scarcity, highlighting its potential for scientific applications where limited data and reliable uncertainty quantification are critical.

Modeling Unknown Stochastic Dynamical System Subject to External Excitation

Yuan Chen — 2026-04-06T07:00:00Z

SIAM Journal on Scientific Computing, Volume 48, Issue 2, Page C216-C239, April 2026.
Abstract. We present a numerical method for learning an unknown nonautonomous stochastic dynamical systems, i.e., stochastic systems subject to time dependent excitation or control signals. Our basic assumption is that the governing equations for the stochastic system are unavailable. However, short bursts of signal-to-state data consisting of certain known excitation signals and their corresponding system responses are available. When a sufficient amount of such signal-to-state data are available, our method is capable of learning the unknown dynamics and producing an accurate predictive model for the stochastic responses of the system subject to arbitrary excitation signals not in the training data. Our method has two key components: (1) a local approximation of the training signal-to-state data to transfer the learning into a parameterized form and (2) a generative model to approximate the underlying unknown stochastic flow map in distribution. After presenting the method in detail, we present a comprehensive set of numerical examples to demonstrate the performance of the proposed method, especially for long-term system predictions.

Multiscale Neural Networks for Approximating Green’s Functions

Wenrui Hao — 2026-04-06T07:00:00Z

SIAM Journal on Scientific Computing, Volume 48, Issue 2, Page C240-C270, April 2026.
Abstract. Neural networks (NNs) have been widely used to solve partial differential equations (PDEs) in the applications of physics, biology, and engineering. One effective approach for solving PDEs with a fixed differential operator is learning Green’s functions. However, Green’s functions are notoriously difficult to learn due to their poor regularity, which typically requires larger NNs and longer training times. In this paper, we address these challenges by leveraging multiscale NNs to learn Green’s functions. Through theoretical analysis using multiscale Barron space methods and experimental validation, we show that the multiscale approach significantly reduces the necessary NN size and accelerates training.

Fully Discrete Analysis of the Galerkin POD Neural Network Approximation with Application to 3D Acoustic Wave Scattering

Jürgen Dölz — 2026-04-08T07:00:00Z

SIAM Journal on Scientific Computing, Volume 48, Issue 2, Page C271-C300, April 2026.
Abstract. In this work, we consider the approximation of parametric maps using the so-called Galerkin POD-NN method. This technique combines the computation of a reduced basis via proper orthogonal decomposition (POD) and artificial neural networks (NNs) for the construction of fast surrogates of said parametric maps. In contrast to the existing literature, which has studied the approximation properties of this kind of architecture on a continuous level, we provide a fully discrete error analysis of this approach. More precisely, our estimates also account for discretization errors during the construction of the NN architecture. We consider the number of reduced basis in the approximation of the solution manifold, truncation in the parameter space, and, most importantly, the number of samples in the computation of the reduced space, together with the effect of the use of NNs in the approximation of the reduced coefficients. Following this error analysis, we provide a priori bounds on the required POD tolerance, the resulting POD ranks, and NN parameters to maintain the order of convergence of quasi Monte Carlo sampling techniques. We conclude this work by showcasing the applicability of this method through a practical industrial application: the sound-soft acoustic scattering problem by a parametrically defined scatterer in three physical dimensions. Reproducibility of computational results. This paper has been awarded the “SIAM Reproducibility Badge: Code and data available” as a recognition that the authors have followed reproducibility principles valued by SISC and the scientific computing community. Code and data that allow readers to reproduce the results in this paper are publicly available on the bonndata fileservers at https://doi.org/10.60507/FK2/OJUHWY [].

Exploring Multiple Timescale Dynamics Using Geometric Singular Perturbation-Informed Neural Networks (GSPINNs)

Hamid Mofidi — 2026-04-10T07:00:00Z

SIAM Journal on Scientific Computing, Volume 48, Issue 2, Page C301-C329, April 2026.
Abstract. Multiple timescale systems have long been a subject of extensive study, with geometric singular perturbation theory (GSPT) emerging as a common tool for analyzing such systems. In this work, we present a comprehensive study of ordinary differential equations in the form of boundary value problems that describe multiple timescale systems, focusing specifically on normally hyperbolic systems. We introduce a novel method, geometric singular perturbation-informed neural networks (GSPINNs), which combines the analytical rigor of GSPT with the modeling power of physics-informed neural networks. Our approach converts the system into a connecting problem, allowing for the discovery of detailed dynamical characteristics that provide deeper insights into the system’s behavior. These solutions are then numerically approximated using GSPINNs. We validated our method, GSPINNs, through various challenging problems where traditional methods fail to capture the intricate dynamics of the system. For moderately nonhyperbolic, multiscale problems that are not extremely stiff, GSPINN remains stable enough to yield accurate solutions without resorting to blowup techniques, demonstrating its robustness in such cases. Our analysis and numerical experiments show that GSPINNs produce accurate and stable solutions for stiff fast-slow boundary value problems across a range of challenging examples, whereas conventional finite difference, shooting, and baseline PINN solvers often struggle. This work opens new opportunities for practical applications in modeling and simulating complex systems with varying timescales.

Leveraging Nested MLMC for Sequential Neural Posterior Estimation with Intractable Likelihoods

Xiliang Yang — 2026-04-13T07:00:00Z

SIAM Journal on Scientific Computing, Volume 48, Issue 2, Page C330-C358, April 2026.
Abstract. There is a growing interest in studying sequential neural posterior estimation (SNPE) techniques due to their advantages for simulation-based models with intractable likelihoods. The methods aim to learn the posterior from adaptively proposed simulations using neural network–based conditional density estimators. As an SNPE technique, the automatic posterior transformation (APT) method proposed by Greenberg, Nonnenmacher, and Macke in 2019 performs well and scales to high-dimensional data. However, the APT method requires computing the expectation of the logarithm of an intractable normalizing constant, i.e., a nested expectation. Although atomic proposals were used to render an analytical normalizing constant, it remains challenging to analyze the convergence of learning. In this paper, we reformulate APT as a nested estimation problem. Building on this, we construct several multilevel Monte Carlo (MLMC) estimators for the loss function and its gradients to accommodate different scenarios, including two unbiased estimators, and a biased estimator that trades a small bias for reduced variance and controlled runtime and memory usage. We also provide convergence results of stochastic gradient descent to quantify the interaction of the bias and variance of the gradient estimator. Numerical experiments for approximating complex posteriors with multimodality in moderate dimensions are provided to examine the effectiveness of the proposed methods.

A Deep Solver for Backward Stochastic Volterra Integral Equations

Kristoffer Andersson — 2026-04-16T07:00:00Z

SIAM Journal on Scientific Computing, Volume 48, Issue 2, Page C359-C384, April 2026.
Abstract. We present the first deep-learning solver for backward stochastic Volterra integral equations (BSVIEs) and their fully coupled forward-backward variants. The method trains a neural network to approximate the two solution fields in a single stage, avoiding the use of nested time-stepping cycles that limit classical algorithms. For the decoupled case, we prove a nonasymptotic error bound composed of an a posteriori residual plus the familiar square root dependence on the time step. Numerical experiments are consistent with this rate and reveal two key properties: scalability, in the sense that accuracy remains stable from low dimension up to 500 spatial variables while GPU batching keeps wall-clock time nearly constant; and generality, since the same method handles coupled systems whose forward dynamics depend on the backward solution. These results open practical access to a family of high-dimensional, time-inconsistent problems in stochastic control and quantitative finance. Reproducibility of computational results. This paper has been awarded the “SIAM Reproducibility Badge: Code and Data Available” as a recognition that the authors have followed reproducibility principles valued by SISC and the scientific computing community. Code and data that allow readers to reproduce the results in this paper are available at https://github.com/AlessandroGnoatto/DeepBSVIE and in the supplementary materials (DeepBSVIE-main.zip [1.76MB]), linked from the main article webpage.

A Paired Autoencoder Framework for Inverse Problems via Bayes Risk Minimization

Emma Hart — 2026-04-17T07:00:00Z

SIAM Journal on Scientific Computing, Volume 48, Issue 2, Page C385-C414, April 2026.
Abstract. In this work, we introduce a data-driven approach for solving inverse problems that exploits technologies from machine learning, in particular autoencoder network structures. We consider a paired autoencoder framework, in which two autoencoders are used to efficiently represent the input and target spaces separately and optimal mappings are learned between latent spaces, thus enabling forward and inverse surrogate mappings. We focus on interpretations using Bayes risk and empirical Bayes risk minimization, and we provide various theoretical results for linear cases, drawing connections to existing works on low-rank matrix approximations. Similar to end-to-end approaches, our paired approach creates a surrogate model for forward propagation and regularized inversion. However, by separately addressing tasks of compression and learning the inverse operator, our approach can outperform existing methods that learn the end-to-end mapping directly. Furthermore, we show that cheaply computable evaluation metrics are available through this framework and can be used to predict whether the solution for a new sample should be predicted well. Reproducibility of computational results. This paper has been awarded the “SIAM Reproducibility Badge: Code and data available” as a recognition that the authors have followed reproducibility principles valued by SISC and the scientific computing community. Code and data that allow readers to reproduce the results in this paper are available at https://github.com/emmahart2000/PAIR and in the supplementary materials (PAIR-main.zip [43.9KB]).

Straggler-Tolerant Stationary Methods for Linear Systems

Vassilis Kalantzis — 2025-01-17T08:00:00Z

SIAM Journal on Scientific Computing, Ahead of Print.
Abstract. In this paper, we consider the iterative solution of sparse systems of linear algebraic equations under the condition that sparse matrix-vector products with the coefficient matrix are computed only partially. At the same time, noncomputed entries are set to zero. We assume that both the number of computed entries and their associated row index set are random variables, with the row index set sampled uniformly given the number of computed entries. This model of computations is prevalent to that realized in hybrid cloud computing architectures following the controller-worker distributed model under the influence of straggling workers. We propose a randomized Richardson iterative scheme and a randomized Chebyshev semi-iterative method within this model and prove the sufficient conditions for their convergence in expectation. Numerical experiments verify the presented theoretical results as well as the effectiveness of the proposed schemes on a few sparse matrix problems.

Asynchronous Semi-iterative Methods and the Asynchronous Chebyshev Method with Multigrid Preconditioning

Jordi Wolfson-Pou — 2025-03-21T07:00:00Z

SIAM Journal on Scientific Computing, Ahead of Print.
Abstract. This paper considers solving linear systems of equations by asynchronous versions of semi-iterative methods which involve one or more parameters. Semi-iterative methods can converge very rapidly when these parameters are chosen well. We first observe that the best parameters for an asynchronous semi-iterative method can be quite different from the best parameters for the corresponding standard, synchronous method. We then present an asynchronous version of the Chebyshev semi-iterative method. As a second-order method, the Chebyshev method is more sensitive to asynchronous execution than first-order semi-iterative methods. This sensitivity can be reduced by different choices of its parameters (e.g., when underestimating the spectrum of the matrix), and also by preconditioning. This motivates a major part of this paper, which is the development of an asynchronous additive multigrid preconditioner for the asynchronous Chebyshev method that is based on iteratively solving an extended, semidefinite system. We demonstrate a shared memory parallel implementation that uses this extended matrix efficiently in matrix-free fashion.

Improving Greedy Algorithms for Rational Approximation

James H. Adler — 2025-04-07T07:00:00Z

SIAM Journal on Scientific Computing, Ahead of Print.
Abstract. When developing robust preconditioners for multiphysics problems, fractional functions of the Laplace operator often arise and need to be inverted. Rational approximation in the uniform norm can be used to convert inverting those fractional operators into inverting a series of shifted Laplace operators. Care must be taken in the approximation so that the shifted Laplace operators remain symmetric positive definite, making them better conditioned. In this work, we study two greedy algorithms for finding rational approximations to such fractional operators. The first algorithm improves the orthogonal greedy algorithm discussed in [Y. Li, L. Zikatanov, and C. Zuo, SIAM J. Sci. Comput., 46 (2024), pp. S68–S87] by adding one minimization step in the uniform norm to the procedure. The second approach employs the weak Chebyshev greedy algorithm in the uniform norm. Both methods yield nonincreasing error. Numerical results confirm the effectiveness of our proposed algorithms, which are also flexible and applicable to other approximation problems. Moreover, with effective rational approximations to the fractional operator, the resulting algorithms show good performance in preconditioning a Darcy–Stokes coupled problem.

An Optimization-Based Coupling of Reduced Order Models with an Efficient Reduced Adjoint Basis Generation Approach

Elizabeth Hawkins — 2025-04-15T07:00:00Z

SIAM Journal on Scientific Computing, Ahead of Print.
Abstract. Optimization-based coupling (OBC) is an attractive alternative to traditional Lagrange multiplier approaches in multiple modeling and simulation contexts. However, application of OBC to time-dependent problems has been hindered by the computational cost of finding the stationary points of the associated Lagrangian, which requires primal and adjoint solves. This issue can be mitigated by using OBC in conjunction with computationally efficient reduced order models (ROMs). To demonstrate the potential of this combination, in this paper, we develop an optimization-based ROM-ROM coupling for a transient advection-diffusion transmission problem. We pursue the “optimize-then-reduce” path toward solving the minimization problem at each time step and solve reduced space adjoint system of equations, where the main challenge in this formulation is the generation of adjoint snapshots and reduced bases for the adjoint systems required by the optimizer. One of the main contributions of the paper is a new technique for an efficient adjoint snapshot collection for gradient-based optimizers in the context of optimization-based ROM-ROM couplings. We present numerical studies demonstrating the accuracy of the approach along with comparison between various approaches for selecting a reduced order basis for the adjoint systems, including decay of snapshot energy, average iteration counts, and timings.

Generalized Optimal AMG Convergence Theory for Nonsymmetric and Indefinite Problems

Ahsan Ali — 2025-05-06T07:00:00Z

SIAM Journal on Scientific Computing, Ahead of Print.
Abstract. Algebraic multigrid (AMG) is known to be an effective solver for many sparse symmetric positive definite (SPD) linear systems. For SPD systems, the convergence theory of AMG is well-understood in terms of the [math]-norm, but in a nonsymmetric setting such an energy norm is nonexistent. For this reason, convergence of AMG for nonsymmetric systems of equations remains an open area of research. A particular aspect missing from theory of nonsymmetric and indefinite AMG is the incorporation of general relaxation schemes. In the SPD setting, the classical form of optimal AMG interpolation provides a useful insight in determining the best possible two-grid convergence rate of a method based on an arbitrary symmetrized relaxation scheme. In this work, we discuss a generalization of the optimal AMG convergence theory targeting nonsymmetric problems, using a certain matrix-induced orthogonality of the left and right eigenvectors of a generalized eigenvalue problem relating the system matrix and relaxation operator. We show that using this generalization of the optimal convergence theory, one can obtain a measure of the spectral radius of the two-grid error transfer operator that is mathematically equivalent to the derivation in the SPD setting for optimal interpolation, which instead uses norms. In addition, this generalization of the optimal AMG convergence theory can be further extended for symmetric indefinite problems, such as those arising from saddle-point systems so that one can obtain a precise convergence rate of the resulting two-grid method based on optimal interpolation. We provide supporting numerical examples of the convergence theory for nonsymmetric advection-diffusion problems, the two-dimensional Dirac equation motivated by [math]-symmetry, and the mixed Darcy flow problem corresponding to a saddle-point system.

Parameter Selection by GCV and a [math] Test within Iterative Methods for [math]-Regularized Inverse Problems

Brian Sweeney — 2025-05-06T07:00:00Z

SIAM Journal on Scientific Computing, Ahead of Print.
Abstract. [math] regularization is used to preserve edges or enforce sparsity in a solution to an inverse problem. We investigate the split Bregman and the majorization-minimization iterative methods that turn this nonsmooth minimization problem into a sequence of steps that include solving an [math]-regularized minimization problem. We consider selecting the regularization parameter in the inner generalized Tikhonov regularization problems that occur at each iteration in these [math] iterative methods. The generalized cross validation method and [math] degrees of freedom test are extended to these inner problems. In particular, for the [math] test this includes extending the [math] result for problems in which the regularization operator has more rows than columns and showing how to use the [math]-weighted generalized inverse to estimate prior information at each inner iteration. Numerical experiments for image deblurring problems demonstrate that it is more effective to select the regularization parameter automatically within the iterative schemes than to keep it fixed for all iterations. Moreover, an appropriate regularization parameter can be estimated in the early iterations and fixed to convergence. Reproducibility of computational results. This paper has been awarded the “SIAM Reproducibility Badge: Code and data available” as a recognition that the authors have followed reproducibility principles valued by SISC and the scientific computing community. Code and data that allow readers to reproduce the results in this paper are available at https://github.com/briansweeney397/ParamEstSBMM and in the supplementary materials (ParamEstSBMM-main.zip [192KB]).

Asymptotic Convergence of Restarted Anderson Acceleration for Certain Normal Linear Systems

Oliver A. Krzysik — 2025-05-07T07:00:00Z

SIAM Journal on Scientific Computing, Ahead of Print.
Abstract. Anderson acceleration (AA) is widely used for accelerating the convergence of an underlying fixed-point iteration [math], [math], with [math], [math]. Despite AA’s widespread use, relatively little is understood theoretically about the extent to which it may accelerate the underlying fixed-point iteration. To this end, we analyze a restarted variant of AA with a restart size of one, a method closely related to GMRES(1). We consider the case of [math] with matrix [math] either symmetric or skew-symmetric. For both classes of [math] we compute the worst-case root-average asymptotic convergence factor of the AA method, partially relying on conjecture in the symmetric setting, proving that it is strictly smaller than that of the underlying fixed-point iteration. For symmetric [math], we show that the AA residual iteration corresponds to a fixed-point iteration for solving an eigenvector-dependent nonlinear eigenvalue problem (NEPv), and we show how this can result in the convergence factor strongly depending on the initial iterate, which we quantify exactly in certain special cases. Conversely, for skew-symmetric [math] we show that the AA residual iteration is closely related to a power iteration for [math], and how this results in the convergence factor being independent of the initial iterate. Supporting numerical results are given, which also indicate the theory is applicable to the more general setting of nonlinear [math] with Jacobian at the fixed point that is symmetric or skew-symmetric. Reproducibility of computational results. This paper has been awarded the “SIAM Reproducibility Badge: Code and data available” as a recognition that the authors have followed reproducibility principles valued by SISC and the scientific computing community. Code and data that allow readers to reproduce the results in this paper are available at https://github.com/okrzysik/rAA1-paper-code and in the supplementary material (M167226_SM.pdf [2.60MB], M167226_SM.zip [8.22MB]).

Inner-Product Free Krylov Methods for Large-Scale Inverse Problems

Ariana N. Brown — 2025-08-01T07:00:00Z

SIAM Journal on Scientific Computing, Ahead of Print.
Abstract. In this study, we introduce two new Krylov subspace methods for solving rectangular large-scale linear inverse problems. The first approach is a modification of the Hessenberg iterative algorithm that is based on an LU factorization and is therefore referred to as the least squares LU (LSLU) method. The second approach incorporates Tikhonov regularization in an efficient manner; we call this the Hybrid LSLU method. Both methods are inner-product free, making them advantageous for high performance computing and mixed-precision arithmetic. Theoretical findings and numerical results show that Hybrid LSLU can be effective in solving large-scale inverse problems and has performance comparable to that of existing iterative projection methods.

Efficient Solution of Fully Implicit Runge–Kutta Methods for Linear Wave Equations

Aman Rani — 2025-08-11T07:00:00Z

SIAM Journal on Scientific Computing, Ahead of Print.
Abstract. We focus on time integration of time-dependent linear hyperbolic PDEs, where implicit Runge–Kutta (IRK) time stepping is an increasingly popular technique. However, this approach gives rise to linear block systems with a Kronecker product structure, involving the Butcher table and the mass and stiffness matrices. Taking the wave equation as the canonical example, the system arising from IRK time stepping is highly ill-conditioned but we can exploit the block structure. If [math] is the number of degrees of freedom for the discretization of the Laplace operator, then the resulting system matrix is a block [math] matrix where each block is of size [math], and [math] is the number of IRK stages. We reformulate the large [math] block structured linear system as a Sylvester matrix equation. This leads to [math] separate systems of order [math], and these smaller systems are efficiently handled with the subsolves replaced by a single AMG V-cycle. We demonstrate the effectiveness of our approach on a two-dimensional wave problem. Our experiments show that our approach not only reduces runtime but also requires fewer AMG V-cycles compared to traditional methods. As the number of Runge–Kutta stages increases and the mesh is refined, the Sylvester approach proves to be at least twice as fast as other existing methods, while also requiring fewer AMG V-cycles. We also introduce a block lower triangular preconditioner based on minimization of [math] over the lower triangular matrices [math] ([math] being the Butcher table), which improves on an existing method based on minimization of [math].

Monolithic Multigrid Preconditioners for High-Order Discretizations of Stokes Equations

Alexey Voronin — 2025-08-13T07:00:00Z

SIAM Journal on Scientific Computing, Ahead of Print.
Abstract. This work introduces and assesses the efficiency of a monolithic [math]MG multigrid framework designed for high-order discretizations of stationary Stokes systems using Taylor–Hood and Scott–Vogelius elements. The proposed approach integrates coarsening in both approximation order ([math]) and mesh resolution ([math]) to address the computational and memory efficiency challenges that are often encountered in conventional high-order numerical simulations. Our numerical results reveal that [math]MG offers significant improvements over traditional spatial-coarsening-only multigrid ([math]MG) techniques for problems discretized with Taylor–Hood elements across a variety of problem sizes and discretization orders. In particular, the [math]MG method exhibits superior performance in reducing setup and solve times, particularly when dealing with higher discretization orders and unstructured problem domains. For Scott–Vogelius discretizations, while monolithic [math]MG delivers low iteration counts and competitive solve phase timings, it exhibits a discernibly slower setup phase when compared to a multilevel (nonmonolithic) full-block-factorization (FBF) preconditioner where [math]MG is employed only for the velocity unknowns. This is primarily due to the setup costs of the larger mixed-field relaxation patches with monolithic [math]MG versus the patch setup costs with a single unknown type for FBF. Reproducibility of computational results. This paper has been awarded the “SIAM Reproducibility Badge: Code and data available” as a recognition that the authors have followed reproducibility principles valued by SISC and the scientific computing community. Code and data that allow readers to reproduce the results in this paper are available at https://github.com/Alexey-Voronin/HighOrderStokesMG and in the supplementary materials (HighOrderStokesMG-main.zip [2.02MB]).

Multilevel Regularized Newton Methods with Fast Convergence Rates

Nick Tsipinakis — 2025-08-13T07:00:00Z

SIAM Journal on Scientific Computing, Ahead of Print.
Abstract. We introduce new multilevel methods for solving large-scale unconstrained optimization problems. Specifically, the philosophy of multilevel methods is applied to Newton-type methods that regularize the Newton subproblem using second-order information from a coarse (low dimensional) subproblem. The new regularized multilevel methods provably converge from any initialization point and enjoy faster convergence rates than gradient descent. In particular, for arbitrary functions with Lipschitz continuous Hessians, we show that their convergence rate interpolates between the rate of gradient descent and that of the cubic Newton method. If, additionally, the objective function is assumed to be convex, then the proposed method converges with the fast [math] rate. Hence, since the updates are generated using a coarse model in low dimensions, the theoretical results of this paper significantly speed up the convergence of Newton-type or preconditioned gradient methods in practical applications. Preliminary numerical results suggest that the proposed multilevel algorithms are significantly faster than current state-of-the-art methods.

Polynomial Approximation to the Inverse of a Large Matrix

Mark Embree — 2025-08-25T07:00:00Z

SIAM Journal on Scientific Computing, Ahead of Print.
Abstract. The inverse of a large matrix can often be accurately approximated by a polynomial of degree significantly lower than the order of the matrix. The iteration polynomial generated by a run of the GMRES algorithm is a good candidate, and its approximation to the inverse often seems to track the accuracy of the GMRES iteration. We investigate the quality of this approximation through theory and experiment, noting the practical need to add copies of some polynomial terms to improve stability. To mitigate storage and orthogonalization costs, other approaches have appeal, such as polynomial preconditioned GMRES and deflation of problematic eigenvalues. Applications of such polynomial approximations include solving systems of linear equations with multiple right-hand sides (where the solutions to subsequent problems come simply by multiplying the polynomial against the new right-hand sides) and variance reduction in multilevel Monte Carlo methods.

Toward an Algebraic Multigrid Method for the Indefinite Helmholtz Equation

Robert D. Falgout — 2025-08-28T07:00:00Z

SIAM Journal on Scientific Computing, Ahead of Print.
Abstract. It is well-known that multigrid methods are very competitive in solving a wide range of SPD problems. However, achieving such performance for non-SPD matrices remains an open problem. In particular, three main issues may arise when solving a Helmholtz problem: some eigenvalues may be negative or even complex, requiring the choice of an adapted smoother for capturing them, and because the near-kernel space is oscillatory, the geometric smoothness assumption cannot be used to build efficient interpolation rules. Moreover, the coarse correction is not equivalent to a projection method since the indefinite matrix does not define a norm. We present some investigations about designing a method that converges in a constant number of iterations with respect to the wavenumber. The method builds on an ideal reduction-based framework and related theory for SPD matrices to improve an initial least squares minimization coarse selection operator formed from a set of smoothed random vectors. A new coarse correction is proposed to minimize the residual in an appropriate norm for indefinite problems. We also present numerical results at the end of the paper.

Error Analysis and Parallel Scaling Study of a Parareal Parallel-in-Time Integration Algorithm for Particle-in-Fourier Schemes

Sriramkrishnan Muralikrishnan — 2025-10-07T07:00:00Z

SIAM Journal on Scientific Computing, Ahead of Print.
Abstract. We propose a parareal based time parallelization scheme in the phase-space for the particle-in-Fourier (PIF) discretization of the Vlasov–Poisson system used in kinetic plasma simulations. We use PIF with a coarse tolerance for the nonuniform fast Fourier transforms, or the standard particle-in-cell scheme, combined with temporal coarsening, as coarse propagators. This is different from the typical spatial coarsening of particles and/or Fourier modes for parareal, which are not possible or effective for PIF schemes. We perform an error analysis of the algorithm and verify the results numerically with Landau damping, two-stream instability, and Penning trap test cases in 3D-3V. We also implement the space-time parallelization of the PIF schemes in the open-source, performance-portable library IPPL and conduct scaling studies up to 1536 A100 GPUs on the JUWELS booster supercomputer. The space-time parallelization utilizing the parareal algorithm for the time parallelization provides up to 4–6 times speedup compared to spatial parallelization alone and achieves a push rate of around 1 billion particles per second for the benchmark plasma mini-apps considered. Reproducibility of computational results. This paper has been awarded the “SIAM Reproducibility Badge: Code and data available” as a recognition that the authors have followed reproducibility principles valued by SISC and the scientific computing community. Code and data that allow readers to reproduce the results in this paper are available at https://github.com/srikrrish/ippl/tree/parapif-paper-v1.0.0 and in the supplementary materials (ippl-parapif-paper-v1.0.0.zip [589KB]).

Parallel-in-Time Solution of Hyperbolic PDE Systems via Characteristic-Variable Block Preconditioning

O. A. Krzysik — 2025-10-08T07:00:00Z

SIAM Journal on Scientific Computing, Ahead of Print.
Abstract. We consider the parallel-in-time solution of both linear and nonlinear hyperbolic partial differential equation (PDE) systems in one spatial dimension. In the nonlinear setting, the discretized equations are solved with a preconditioned residual iteration based on a global linearization. The linear(ized) equation systems are approximately solved parallel-in-time using a block preconditioner applied in the characteristic variables of the underlying linear(ized) hyperbolic PDE. This change of variables is motivated by the observation that intervariable coupling between characteristic variables is weak, at least locally where spatio-temporal variations in the eigenvectors of the associated flux Jacobian are sufficiently small, while that between the original variables is not. For an [math]-dimensional system of PDEs, applying the preconditioner consists of solving a sequence of [math] scalar linear(ized)-advection-like problems, each associated with a different characteristic wave-speed in the underlying linear(ized) PDE. We approximately solve these linear advection problems using multigrid reduction-in-time (MGRIT); however, any other suitable parallel-in-time method could be used. Numerical examples are shown for the (linear) acoustics equations in heterogeneous media and for the (nonlinear) shallow water equations and Euler equations of gas dynamics with shocks and rarefactions. For many test problems, the solver converges in just a handful of iterations and with mesh-independent convergence rates. Reproducibility of computational results. This paper has been awarded the “SIAM Reproducibility Badge: Code and data available” as a recognition that the authors have followed reproducibility principles valued by SISC and the scientific computing community. Code and data that allow readers to reproduce the results in this paper are available at https://github.com/okrzysik/pit-nonlinear-hyperbolic and in the supplementary materials (M167331_SM.pdf [1.52MB], pit-nonlinear-hyperbolic-master.zip [35.6MB]).

An Adaptive Newton-Based Free-Boundary Grad–Shafranov Solver

Daniel A. Serino — 2025-10-09T07:00:00Z

SIAM Journal on Scientific Computing, Ahead of Print.
Abstract. Equilibria in magnetic confinement devices result from force balancing between the Lorentz force and the plasma pressure gradient. In an axisymmetric configuration like a tokamak, such an equilibrium is described by an elliptic equation for the poloidal magnetic flux, commonly known as the Grad–Shafranov equation. It is challenging to develop a scalable and accurate free-boundary Grad–Shafranov solver, since it is a fully nonlinear optimization problem that simultaneously solves for the magnetic field coil current outside the plasma to control the plasma shape. In this work, we develop a Newton-based free-boundary Grad–Shafranov solver using adaptive finite elements and preconditioning strategies. The free-boundary interaction leads to the evaluation of a domain-dependent nonlinear form of which its contribution to the Jacobian matrix is achieved through shape calculus. The optimization problem aims to minimize the distance between the plasma boundary and specified control points while satisfying two nontrivial constraints, which correspond to the nonlinear finite element discretization of the Grad–Shafranov equation and a constraint on the total plasma current involving a nonlocal coupling term. The linear system is solved by a block factorization, and AMG is called for subblock elliptic operators. The unique contributions of this work include the treatment of a global constraint, preconditioning strategies, nonlocal reformulation, and the implementation of adaptive finite elements. It is found that the resulting Newton solver is robust, successfully reducing the nonlinear residual to 1e-6 and lower in a small handful of iterations while addressing the challenging case to find a Taylor state equilibrium where conventional Picard-based solvers fail to converge. Reproducibility of computational results.This paper has been awarded the “SIAM Reproducibility Badge: Code and data available” as a recognition that the authors have followed reproducibility principles valued by SISC and the scientific computing community. Code and data that allow readers to reproduce the results in this paper are available at https://github.com/mfem/mfem/tree/tds-gs and in the supplementary materials (mfem-tds-gs.zip [18.0MB]).

Transformed Primal-Dual Methods with Variable Preconditioners

Long Chen — 2025-10-13T07:00:00Z

SIAM Journal on Scientific Computing, Ahead of Print.
Abstract. This paper introduces a novel transformed primal-dual with variable preconditioners/metrics (TPDv) algorithm, designed to efficiently solve affine-constrained optimization problems. These problems commonly exist in nonlinear partial differential equations. The TPDv is derived by applying a specially designed transformation matrix to the flow of the primal-dual algorithm, such that the transformed system obtains the strong Lyapunov property. Distinguished from traditional methods, the TPDv iteratively updates time-evolving preconditioning operators, enhancing adaptability. The exponential decay of the TPDv flow and the global linear convergence rates of the TPDv algorithm are established under mild assumptions, where the dependence of the rates on the condition numbers is explicitly specified. Numerical experiments on the Darcy–Forchheimer model and a nonlinear electromagnetic problem demonstrate the algorithm’s superiority over existing methods in terms of iteration numbers and computational time. Reproducibility of computational results. This paper has been awarded the “SIAM Reproducibility Badge: Code and data available” as a recognition that the authors have followed reproducibility principles valued by SISC and the scientific computing community. Code and data that allow readers to reproduce the results in this paper are available at https://github.com/JingrongWei/TPDv_code and in the supplementary materials (TPDv_code-main.zip [40.7KB]).

Efficient Shallow Ritz Method for One-Dimensional Diffusion-Reaction Problems

Zhiqiang Cai — 2025-10-14T07:00:00Z

SIAM Journal on Scientific Computing, Ahead of Print.
Abstract. This paper studies the shallow Ritz method for solving one-dimensional diffusion-reaction problems. The method is capable of improving the order of approximation for nonsmooth problems. By following an approach similar to the one presented in [Z. Cai et al., Efficient Shallow Ritz Method for 1D Diffusion Problems, https://arxiv.org/abs/2404.17750, 2024], we present a damped block Newton (dBN) method to achieve nearly optimal order of approximation. The dBN method optimizes the Ritz functional by alternating between the linear and nonlinear parameters of the shallow ReLU neural network. For diffusion-reaction problems, new difficulties arise: (1) for the linear parameters, the mass matrix is dense and even more ill-conditioned than the stiffness matrix, and (2) for the nonlinear parameters, the Hessian matrix is dense and may be singular. This paper addresses these challenges, resulting in a dBN method with computational cost of [math]. The ideas presented for diffusion-reaction problems can also be applied to least-squares approximation problems. For both applications, starting with the nonlinear parameters as a uniform partition, numerical experiments show that the dBN method moves the mesh points to nearly optimal locations.

Convergence Analysis of the Alternating Anderson–Picard Method for Nonlinear Fixed-Point Problems

Xue Feng — 2025-10-28T07:00:00Z

SIAM Journal on Scientific Computing, Ahead of Print.
Abstract. Anderson acceleration (AA) has been widely used to solve nonlinear fixed-point problems due to its rapid convergence. This work focuses on a variant of AA in which multiple Picard iterations are performed between each AA step, referred to as the Alternating Anderson–Picard (AAP) method. Despite introducing more “slow” Picard iterations, this method has been shown to be efficient and even more robust in both linear and nonlinear cases. However, there is a lack of theoretical analysis for AAP in the nonlinear case. In this paper, we address this gap by establishing the equivalence between AAP and a multisecant-GMRES method that uses GMRES to solve a multisecant linear system at each iteration. From this perspective, we show that AAP “converges” to the Newton-GMRES method. Specifically, as the residual approaches zero, the multisecant matrix, the approximate Jacobian inverse, the search direction, and the optimization gain of AAP converge to their counterparts in the Newton-GMRES method. These connections provide insights for analyzing the asymptotic convergence properties of AAP. Consequently, we show that AAP is locally [math]-linear convergent and provide an upper bound for the convergence factor of AAP. To validate the theoretical results, numerical examples are provided. Reproducibility of computational results.This paper has been awarded the “SIAM Reproducibility Badge: Code and data available” as a recognition that the authors have followed reproducibility principles valued by SISC and the scientific computing community. Code and data that allow readers to reproduce the results in this paper are available at https://github.com/xue1993/AAP.git and in the supplementary materials (AAP-main.zip [2.06MB]).

On a Randomized Small-Block Lanczos Method for Large-Scale Null Space Computations

Daniel Kressner — 2026-04-02T07:00:00Z

SIAM Journal on Scientific Computing, Ahead of Print.
Abstract. Computing the null space of a large sparse matrix [math] is a challenging computational problem, especially if the nullity—the dimension of the null space—is not small. When applying a block Lanczos method to [math] for this purpose, conventional wisdom suggests to use a block size [math] that is not smaller than the nullity. In this work, we show how randomness can be utilized to allow for smaller [math] without sacrificing convergence or reliability. Even [math], corresponding to the standard single-vector Lanczos method, becomes a safe choice. This is achieved by using a small random diagonal perturbation, which moves the zero eigenvalues of [math] away from each other, and a random initial guess. We analyze the effect of the perturbation on the attainable quality of the null space and derive convergence results that establish robust convergence for [math]. As demonstrated by our numerical experiments, a smaller block size combined with restarting and partial reorthogonalization results in reduced memory requirements and computational effort. It also allows for the incremental computation of the null space without requiring a priori knowledge of the nullity. Our algorithm is best suited for situations when the nullity of [math] is moderate. Reproducibility of computational results. This paper has been awarded the “SIAM Reproducibility Badge: Code and data available” as a recognition that the authors have followed reproducibility principles valued by SISC and the scientific computing community. Code and data that allow readers to reproduce the results in this paper are available at https://github.com/nShao678/Nullspace-code and in the supplementary materials (Nullspace-code-main.zip [10.7MB]).

Distributed Memory Parallel Adaptive Tensor-Train Cross Approximation

Tianyi Shi — 2026-04-17T07:00:00Z

SIAM Journal on Scientific Computing, Ahead of Print.
Abstract. The tensor-train (TT) format is a data-sparse tensor representation commonly used in high dimensional function approximations arising from computational and data sciences. Various sequential and parallel TT decomposition algorithms have been proposed for different tensor inputs and assumptions. In this paper, we propose subtensor parallel adaptive TT cross, which partitions a tensor onto distributed memory machines with multidimensional process grids and constructs an TT approximation iteratively with tensor elements. We derive two iterative formulations for pivot selection and TT core construction under the distributed memory setting, conduct communication and scaling analysis of the algorithm, and illustrate its performance with multiple test experiments. These include up to 6D Hilbert tensors and tensors constructed from Maxwellian distribution functions that arise in kinetic theory. Our results demonstrate significant accuracy with greatly reduced storage requirements via the TT cross approximation. Furthermore, we demonstrate good to optimal strong and weak scaling performance for the proposed parallel algorithm. Reproducibility of computational results. This paper has been awarded the “SIAM Reproducibility Badge: Code and data available” as a recognition that the authors have followed reproducibility principles valued by SISC and the scientific computing community. Code and data that allow readers to reproduce the results in this paper are available at https://github.com/dhayes95/Cross and in the supplementary materials (Cross-main.zip [15KB]).