Abstract. Differential equations (DEs) are widely used to model the evolution of natural or engineered systems. Traditionally, the parameters in DEs are adjusted to fit data obtained through system observations. However, in fields such as politics, economics, and biology, the observed data points are often collected independently (i.e., repeated cross-sectional (RCS) data). In this paper, we first demonstrate that conventional optimization techniques struggle to accurately estimate DE parameters when RCS data exhibit various heterogeneities, leading to a significant loss of information. To address this issue, we propose a new estimation method, the emulator-informed deep-generative model (EIDGM), designed to handle RCS data. Specifically, EIDGM integrates a physics-informed neural network–based emulator that immediately generates DE solutions and a Wasserstein-generative adversarial network–based parameter generator that can effectively mimic the RCS data. We then evaluated the effectiveness of EIDGM across various models, including exponential and logistic growth models, and the Lorenz system, demonstrating its superior ability to accurately capture parameter distributions. In addition, we applied EIDGM to real-world datasets, successfully capturing diverse shapes of parameter distributions. This result highlights that EIDGM can be applied to model a wide range of systems with limited data availability, and our approach can be significantly extended to uncover the operating principles of systems based on limited data. 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/CHWmath/EIDGM/tree/main and in the supplementary materials ( EIDGM-main.zip [19.5MB]). [Formula: see text]
Abstract. Shock waves in gas dynamics feature jump discontinuities that hinder numerical simulations. Viscous regularizations are prone to excessive dissipation of fine-scale structures. In this work, we propose the first inviscid regularization of the multidimensional Euler equation based on ideas from semidefinite programming, information geometry, geometric hydrodynamics, and nonlinear elasticity. The Lagrangian flow maps of Euler solutions are a dynamical system on the manifold of diffeomorphisms. We observe that shock formation arises from the manifold’s geodesic incompleteness. Our regularization embeds it into an ambient space equipped with the information geometry of the logarithmic barrier function. Thus, the diffeomorphism manifold inherits a geodesically complete geometry. The resulting regularized conservation law replaces shocks with smooth profiles without affecting oscillatory structures. One- and two-dimensional numerical experiments show its practical potential to enable higher-order methods without explicit shock capturing. While we focus on the barotropic Euler equations for concreteness and simplicity of exposition, our regularization easily extends to more general Euler and Navier–Stokes-type equations. Our approach regularizes the Wasserstein geometry of the mass density with its information geometry. The former captures the natural trajectories of physical particles and the latter that of statistical estimators. Information geometric regularization accounts for the mass density’s dual nature as a statistical/computational tool summarizing the motion of physical particles. Thus, our work is a starting point for information geometric mechanics that views solutions of continuum mechanical PDEs as parameters of statistical models for unresolved scales and uses their information geometry to evolve them in 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/f-t-s/information_geometric_regularization_of_barotropic_euler and in the supplementary materials ( information_geometric_regularization_of_barotropic_euler-main.zip [86.2KB]). [Formula: see text]
Abstract. A wide range of applications necessitates solving large-scale ill-posed problems contaminated by noise, which are commonly formulated as linear approximation problems (possibly overdetermined or underdetermined). Krylov subspace regularization methods are particularly advantageous in this context, as they rely solely on matrix-vector multiplication. Among the most widely used techniques are least squares QR (LSQR) and conjugate gradient least squares (CGLS), both of which can be extended with flexible preconditioning to enforce solution properties such as nonnegativity or sparsity. Flexible LSQR (FLSQR) can also be combined with direct methods to create efficient hybrid approaches. The flexible Golub–Kahan process underlying FLSQR requires two long-term recurrences. Here, we introduce a novel fast flexible Golub–Kahan method that employs one long-term and one short-term recurrence. Using this, we develop the fast flexible LSQR (FaFLSQR) algorithm, which offers comparable computational cost to flexible CGLS (FCGLS) while also supporting hybrid regularization like FLSQR. We analyze the properties of FaFLSQR and prove its mathematical equivalence to FCGLS. Numerical experiments demonstrate that in floating-point arithmetic, FaFLSQR outperforms both FCGLS and FLSQR in terms of computational efficiency.
Abstract. The impressive performance of deep neural networks (DNNs) on a variety of learning tasks has spurred much investigation into improving their training and characterizing the trained functions. Recent work has shown the equivalence of DNNs in the infinite width limit and kernel machines relying on the neural tangent kernel (NTK) at initialization. These results suggest, and experimental evidence corroborates, that kernel machines relying on empirical kernels extracted from trained DNNs can act as surrogates for trained finite-width DNNs. The high computational cost of assembling the NTK, however, makes this approach infeasible in practice. In the current work, we study the performance of the conjugate kernel (CK), an efficient approximation to the NTK. For smooth function and logistic regression, we show that the CK performance is only marginally worse than that of the NTK and, in certain cases, much more superior. In particular, we establish bounds for the test losses, verify them with numerical tests, and identify the regularity of the kernel as the key determinant of performance. We also determine regimes where both kernel machines trained on features extracted from an underlying DNN are demonstrably superior to the latter and use this to suggest a recipe for accelerating DNN performance inexpensively. We present a demonstration of this on foundation models by comparing their performance on a classification task using a conventional technique and our prescription. We also show how our approach can be used to improve physics-informed operator network training as well as convolutional neural network training for vision classification tasks.
Abstract. We propose a particle-based workflow for approximating the time-dependent law of finite-volume discretizations of the Dean–Kawasaki model. After discretization, the state is a nonnegative vector whose total mass is conserved by the finite-volume scheme, and it is therefore supported on a probability simplex. To enable tensor-network density estimation, we map the simplex to an unconstrained Euclidean space using a centered logarithmic transform and then apply a wavelet transform that organizes degrees of freedom by spatial scale. On the transformed variables, we fit the probability density with a functional hierarchical tensor over a wavelet basis, i.e., a hierarchical-Tucker/tree-tensor-network representation of the coefficient tensor of a fixed univariate basis expansion. We illustrate the method on 1D and 2D examples with [Formula: see text] degrees of freedom, including cases with external potentials and pairwise interactions. The method accurately captures the site-wise correlations and other observables of the true model. 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/Xun-Tang123/FHT_for_deans_equation and in the supplementary materials ( FHT_for_deans_equation-main.zip [27.5MB]). [Formula: see text]
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 Nunberg (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.
We introduce a new Approximate Cholesky factorization algorithm and software for solving linear equations in symmetric diagonally dominant matrices with nonpositive off-diagonal entries (SDDM matrices), including Laplacian matrices. We use preconditioned conjugate gradient to solve these systems. Our preconditioner is a variant of the Approximate Cholesky factorization of Kyng and Sachdeva (FOCS 2016). Our factorization approach is simple: we eliminate matrix rows/columns one at a time and use sampling to update the remaining entries of the matrix, approximating complete Cholesky factorization. Unlike Kyng-Sachdeva, we do not prove theoretical guarantees about our algorithm, but we observe experimentally that our algorithm runs much faster. Our algorithm broadly resembles a common approach of using Incomplete Cholesky factorization as a preconditioner in conjugate gradient (ICCG). Unlike earlier approaches such as Kyng-Sachdeva and Incomplete Cholesky variants, our sampled entries always maintain a connected support graph on the neighbors of the eliminated variable. We implement our solver algorithm in Julia and experimentally evaluate its performance when using 1 or 2 samples for each original entry. We compare our single-threaded performance to that of state-of-the-art solvers. We evaluate performance across many problem classes, including all large SDDM matrices from the SuiteSparse collection, and a broad array of programmatically generated instances, with tests ranging up to 200 million nonzeros per system of linear equations. We hope that others will adopt this suite of tests as a benchmark, which we refer to as SDDM2023. Our experiments suggest that our solvers are much more robust than existing solvers for SDDM linear equations, while retaining good performance across all instances. Our code is public, and we detail precisely the tests we run and provide a tutorial on how to replicate them.
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.
Solving multiscale diffusion problems is often computationally expensive due to the spatial and temporal discretization challenges arising from high-contrast coefficients. To address this issue, a partially explicit temporal splitting scheme is proposed. By appropriately constructing multiscale spaces, the spatial multiscale property is effectively captured, and it has been demonstrated that the temporal step size is independent of the contrast. To enhance simulation speed, we propose a parallel algorithm for the multiscale flow problem that leverages the partially explicit temporal splitting scheme. The idea is first to evolve the partially explicit system using a coarse time step size and then correct the solution on each coarse time interval with a fine propagator, for which we consider the all-at-once solver. This procedure is then performed iteratively until convergence. We analyze the stability and convergence of the proposed algorithm. The numerical experiments demonstrate that the proposed algorithm achieves high numerical accuracy for high-contrast problems and converges in a relatively small number of iterations. The number of iterations stays stable as the number of coarse intervals increases, thus improving computational efficiency through parallel processing.
This paper addresses the inverse scattering problem of reconstructing multilayered elastic periodic media using only one-sided near-field measurements. A novel Lippmann--Schwingertype integral equation is established by introducing a two-layer quasiperiodic Green's tensor for the Navier equation, providing a foundation for the well-posedness analysis of the direct problem. Building on this formulation, we propose an innovative inversion algorithm that reconstructs the whole periodic structure using one-sided measurements generated by only a small number of compressional waves. The method combines a Newton-type iterative scheme with total variation (TV) regularization for sharp interface recovery, and incorporates a segmentation-based postprocessing step to improve the accuracy and enforce piecewise constancy. Numerical experiments show that the algorithm accurately reconstructs complex multilayer structures with sharp corners, which are beyond the reach of existing methods even in the simpler Helmholtz case. Notably, even when the measured data contain little to no evanescent wave information, leading to a highly ill-posed inverse problem, experiments demonstrate that our algorithm still provides accurate reconstructions.
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 [6], 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.
Structure-preservation, accuracy, and robustness are critical aims for magnetohydro dynamic (MHD) simulations. In this paper, we propose a high-order, magnetic helicity conservmethod 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 structurepreserving 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.
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.
High-resolution reservoir simulation of multiphase and multicomponent porous media flows, arising from geological carbon sequestration, requires solving large, sparse nonlinear and linear systems efficiently on parallel computers. This multiphysical process involves a compositional flow model with multiple wells and corner-point grid geometry, which demands efficient and scalable solvers optimized for high-performance computing. In the study, we present a nonstandard restricted additive Schwarz preconditioner with Schur complement, called SchurRAS, to mitigate slow convergence and stagnation in solving sparse linear systems. Our hybrid preconditioning approach combines overlapping domain decomposition for parallel scalability with localized Schur complement approximations to enhance solver efficiency on distributed-memory systems equipped with multicore nodes. Moreover, a minimum-type semismooth Newton algorithm is proposed to enhance nonlinear convergence and enforce physical constraints. Large-scale simulations are conducted on geological reservoir benchmarks and real-field carbon sequestration projects. The parallel performance evaluation on a supercomputer demonstrates the scalability of our approach to thousands of processor cores and shows that SchurRAS outperforms standard RAS in both iteration count and computational time.
Neural network tearing and interconnecting methods are proposed and tested for second-order elliptic problems to address the longer training time in neural network solutions for PDEs. The methods are developed for nonoverlapping sub domain partitions of the problem domain, and such nonoveralapping partitions are useful for modeling coefficient discontinuity, multi-physics, and interface discontinuity inside the problem domain. A gradient-based iterative algorithm is first derived for the solution on the sub domain interface, where, at each iteration, local neural network solutions are trained for the provided interface solution values independently, allowing parallel computation, and the interface solution values are updated by using the trained local neural network solutions. The convergence of the iterative algorithm is shown by using well-established finite element tearing and interconnecting methods. The gradient-based iterative algorithm is then further enhanced by proposing preconditioning schemes on the gradient value. The performance of the proposed preconditioning schemes is demonstrated for various test examples.
Abstract. Gaussian process (GP)–based methods for solving PDEs demonstrate great promise by bridging the gap between the theoretical rigor of traditional numerical algorithms and the flexible design of machine learning solvers. The main bottleneck of GP methods lies in the inversion of a covariance matrix, whose cost grows cubically concerning the size of samples. Drawing inspiration from neural networks, we propose a minibatch algorithm combined with GPs to solve nonlinear PDEs. A naive deployment of a stochastic gradient descent method for solving PDEs with GPs is challenging, as the objective function in the requisite minimization problem cannot be depicted as the expectation of a finite-dimensional random function. To address this issue, we employ a minibatch method to the corresponding infinite-dimensional minimization problem over function spaces. The algorithm takes a minibatch of samples at each step to update the GP model. Thus, the computational cost is allotted to each iteration. Using stability analysis and convexity arguments, we show that, in expectation, the iterates of the minibatch algorithm approach a nearly stationary point of the optimal recovery problem used to solve the PDEs with rate [Formula: see text], where [Formula: see text] is the number of iterations and [Formula: see text] is the batch size.
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 eigenmo des. 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 N (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.
Abstract. Spectral methods employing nonstandard polynomial bases, such as Müntz polynomials, have proven effective for accurately solving singular problems. However, due to the absence of theoretical guidance, the key parameter [Formula: see text] that controls the exponents of Müntz polynomials is usually determined empirically through extensive numerical experiments, leading to a time-consuming tuning process. To address this bottleneck, we propose a machine-learning–based optimization framework for the Müntz spectral method. As an illustrative example, we optimize the parameter selection for solving time-fractional partial differential equations (PDEs). Specifically, an artificial neural network (ANN) is employed to predict optimal [Formula: see text] based solely on the time-fractional order [Formula: see text]. The ANN is trained by minimizing the mean solution error on a set of sampled one-dimensional time-fractional convection-diffusion equations. Numerical results for time-fractional PDEs demonstrate that the ANN-based parameter prediction significantly improves the accuracy of the Müntz spectral method. Moreover, the trained ANN generalizes effectively from one-dimensional to two-dimensional cases, highlighting its robustness across spatial dimensions.
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.
