2025-06-19 | | Total: 25
We consider the problem of estimating the spectral norm of a matrix using only matrix-vector products. We propose a new Counterbalance estimator that provides upper bounds on the norm and derive probabilistic guarantees on its underestimation. Compared to standard approaches such as the power method, the proposed estimator produces significantly tighter upper bounds in both synthetic and real-world settings. Our method is especially effective for matrices with fast-decaying spectra, such as those arising in deep learning and inverse problems.
The $h$-version of the finite-element method ($h$-FEM) applied to the high-frequency Helmholtz equation has been a classic topic in numerical analysis since the 1990s. It is now rigorously understood that (using piecewise polynomials of degree $p$ on a mesh of a maximal width $h$) the conditions "$(hk)^p \rho$ sufficiently small" and "$(hk)^{2p} \rho$ sufficiently small" guarantee, respectively, $k$-uniform quasioptimality (QO) and bounded relative error (BRE), where $\rho$ is the norm of the solution operator with $\rho\sim k$ for non-trapping problems. Empirically, these conditions are observed to be optimal in the context of $h$-FEM with a uniform mesh. This paper demonstrates that QO and BRE can be achieved using certain non-uniform meshes that violate the conditions above on $h$ and involve coarser meshes away from trapping and in the perfectly matched layer (PML). The main theorem details how varying the meshwidth in one region affects errors both in that region and elsewhere. One notable consequence is that, for any scattering problem (trapping or nontrapping), in the PML one only needs $hk$ to be sufficiently small; i.e. there is no pollution in the PML. The motivating idea for the analysis is that the Helmholtz data-to-solution map behaves differently depending on the locations of both the measurement and data, in particular, on the properties of billiards trajectories (i.e. rays) through these sets. Because of this, it is natural that the approximation requirements for finite-element spaces in a subset should depend on the properties of billiard rays through that set. Inserting this behaviour into the latest duality arguments for the FEM applied to the high-frequency Helmholtz equation allows us to retain detailed information about the influence of $\textit{both}$ the mesh structure $\textit{and}$ the behaviour of the true solution on local errors in FEM.
This paper delves into the well-posedness and the numerical approximation of non-autonomous stochastic differential algebraic equations (SDAEs) with nonlinear local Lipschitz coefficients that satisfy the more general monotonicity condition called Khasminskii condition. The key challenge is the presence of a singular matrix which makes the numerical integration hard and heavy. To address this challenge, we propose a novel numerical scheme based on semi-implicit method for the drift component of the SDAEs. More precisely we split the drift term as the sum of a linear term and a nonlinear term. The linear part is approximated implicitly, while the nonlinear part is approximated explicitly. The linear component's role is to handle the singularity issues during the numerical integration without the resolution of nonlinear algebraic equations in the constraint equations. This novel scheme is therefore very efficient for SDAEs in high dimension that come after the spatial discretisation of stochastic partial differential algebraic equations (SPDAEs). To prove the pathwise convergence of our novel scheme, we first derive a equivalent scheme called dual scheme, suitable for mathematical analysis and linked to the inherent stochastic differential equation resulting from the elimination of constraints in the initial SDAEs. We prove that our novel scheme converges to the exact solution with rate $\frac{1}{2}-\epsilon$, for arbitrary $\epsilon>0$ in the pathwise sense. Numerical simulations are performed to demonstrate the efficiency of the scheme in high dimension and to show that our theoretical results are in agreement with numerical experiments.
We examine the intrinsic (within the attention head) and extrinsic (amongst the attention heads) structure of the self-attention mechanism in transformers. Theoretical evidence for invariance of the self-attention mechanism to softmax activation is obtained by appealing to paradifferential calculus, (and is supported by computational examples), which relies on the intrinsic organization of the attention heads. Furthermore, we use an existing methodology for hierarchical organization of tensors to examine network structure by constructing hierarchal partition trees with respect to the query, key, and head axes of network 3-tensors. Such an organization is consequential since it allows one to profitably execute common signal processing tasks on a geometry where the organized network 3-tensors exhibit regularity. We exemplify this qualitatively, by visualizing the hierarchical organization of the tree comprised of attention heads and the diffusion map embeddings, and quantitatively by investigating network sparsity with the expansion coefficients of individual attention heads and the entire network with respect to the bi and tri-haar bases (respectively) on the space of queries, keys, and heads of the network. To showcase the utility of our theoretical and methodological findings, we provide computational examples using vision and language transformers. The ramifications of these findings are two-fold: (1) a subsequent step in interpretability analysis is theoretically admitted, and can be exploited empirically for downstream interpretability tasks (2) one can use the network 3-tensor organization for empirical network applications such as model pruning (by virtue of network sparsity) and network architecture comparison.
Recent advances in deep learning makes solving parabolic partial differential equations (PDEs) in high dimensional spaces possible via forward-backward stochastic differential equation (FBSDE) formulations. The implementation of most existing methods requires simulating multiple trajectories of stochastic processes with a small step size of time discretization to ensure accuracy, hence having limited performance, especially when solving on a large time interval. To address such issue, we propose a deep "shotgun method" that does not exploit full trajectories, but only utilizes the data distribution of them. Numerical results including examples with dimensionality up to 10000 demonstrate the competitiveness of the proposed shotgun method in both performance and accuracy.
We study the problem of estimating the diagonal of an implicitly given matrix $\Ab$. For such a matrix we have access to an oracle that allows us to evaluate the matrix quadratic form $ \ub^\top \Ab \ub$. Based on this query oracle, we propose a stochastic diagonal estimation method with random variable $\ub$ drawn from the standard Gaussian distribution. We provide the element-wise and norm-wise sample complexities of the proposed method. Our numerical experiments on different types and dimensions matrices demonstrate the effectiveness of our method and validate the tightness of theoretical results.
In the fields of control theory and machine learning, the dynamic low-rank approximation for large-scale matrices has received substantial attention. Considering the large-scale semilinear stiff matrix differential equations, we propose a dynamic numerical integrator for obtaining low-rank approximations of solutions. We first decompose the differential equation into a stiff linear component and a nonstiff nonlinear term, then employ an exponential integrator along with a dynamic low-rank approach to resolve these subsystems, respectively. Furthermore, the proposed framework naturally extends to rank-adaptation scenarios. Through rigorous validation on canonical stiff matrix differential problems, including spatially discretized Allen-Cahn equations and differential Riccati equations, we demonstrate that the method achieves the theoretically predicted convergence orders. Numerical evidence confirms the robustness and accuracy of the proposed methods.
Hamiltonian particle-based simulations of plasma dynamics are inherently computationally intensive, primarily due to the large number of particles required to obtain accurate solutions. This challenge becomes even more acute in many-query contexts, where numerous simulations must be conducted across a range of time and parameter values. Consequently, it is essential to construct reduced order models from such discretizations to significantly lower computational costs while ensuring validity across the specified time and parameter domains. Preserving the Hamiltonian structure in these reduced models is also crucial, as it helps maintain long-term stability. In this paper, we introduce a nonlinear, non-intrusive, data-driven model order reduction method for the 1D-1V Vlasov--Poisson system, discretized using a Hamiltonian Particle-In-Cell scheme. Our approach relies on a two-step projection framework: an initial linear projection based on the Proper Symplectic Decomposition, followed by a nonlinear projection learned via an autoencoder neural network. The reduced dynamics are then modeled using a Hamiltonian neural network. The offline phase of the method is split into two stages: first, constructing the linear projection using full-order model snapshots; second, jointly training the autoencoder and the Hamiltonian neural network to simultaneously learn the encoder-decoder mappings and the reduced dynamics. We validate the proposed method on several benchmarks, including Landau damping and two-stream instability. The results show that our method has better reduction properties than standard linear Hamiltonian reduction methods.
An efficient procedure using a novel semi-analytical forward solver for identifying heterogeneous and anisotropic elastic parameters from only one full-field measurement is proposed and explored. We formulate the inverse problem as an special energy functional minimization with total variation(TV) regularization. The minimization problem is solved by Adam algorithm, which only requires solving one forward problem and no adjoint problem in each iteration. In order to deal with the irregularity of the elastic regions, the anisotropy and heterogeneity of parameters and potential singularities in forward-modeled issues, a novel semi-analytical forward solver named the direct method of lines is proposed, which discretizes angular variable while preserving analytical solutions along remaining coordinates. To validate the efficacy of our procedure, a series of numerical experiments are implemented subsequently, achieving reliable performance in both forward modeling and the six elastic arguments reconstruction scenarios.
A Fourier transform method is introduced for a class of hybrid time-frequency methods that solve the acoustic scattering problem in regimes where the solution exhibits both highly oscillatory behavior and slow decay in time. This extends the applicability of hybrid time-frequency schemes to domains with trapping regions. A fast sinc transform technique for managing highly oscillatory behavior and long time horizons is combined with a contour integration scheme that improves smoothness properties in the integrand.
To numerically solve the two-dimensional advection equation, we propose a family of fourth- and higher-order semi-Lagrangian finite volume (SLFV) methods that feature (1) fourth-, sixth-, and eighth-order convergence rates, (2) applicability to both regular and irregular domains with arbitrarily complex topology and geometry, (3) ease of handling both zero and nonzero source terms, and (4) the same algorithmic steps for both periodic and incoming penetration conditions. Test results confirm the analysis and demonstrate the accuracy, flexibility, robustness, and excellent conditioning of the proposed SLFV method.
We propose a new nonconforming \(P_1\) finite element method for elliptic interface problems. The method is constructed on a locally anisotropic mixed mesh, which is generated by fitting the interface through a simple connection of intersection points on an interface-unfitted background mesh, as introduced in \cite{Hu2021optimal}. We first establish interpolation error estimates on quadrilateral elements satisfying the regular decomposition property (RDP). Building on this, the main contribution of this work is a novel consistency error analysis for nonconforming elements, which removes the quasi-regularity assumption commonly required in existing approaches. Numerical results confirm the theoretical convergence rates and demonstrate the robustness and accuracy of the proposed method.
There exist elegant methods of aligning point clouds in $\mathbb R^3$. Unfortunately, these methods rely on the positive definite property of the Euclidean metric, and do not easily extend to the indefinite Minkowski metric. In this paper, we propose two solutions to the following problem: given inertial reference frames $A$ and $B$, and given (possibly noisy) measurements of a set of 4-vectors $\{v_i\}$ made in those reference frames with components $\{v_{A,i}\}$ and $\{v_{B,i}\}$, find the optimal Lorentz transformation $\Lambda$ such that $\Lambda v_{A,i}=v_{B,i}$. The method we outline is conceptually simple and easily extends to alignment problems in other matrix Lie groups.
Meshfree methods, including the reproducing kernel particle method (RKPM), have been widely used within the computational mechanics community to model physical phenomena in materials undergoing large deformations or extreme topology changes. RKPM shape functions and their derivatives cannot be accurately integrated with the Gauss-quadrature methods widely employed for the finite element method (FEM) and typically require sophisticated nodal integration techniques, preventing them from easily being implemented in existing FEM software. Interpolation-based methods have been developed to address similar problems with isogeometric and immersed boundary methods, allowing these techniques to be implemented within open-source finite element software. With interpolation-based methods, background basis functions are represented as linear combinations of Lagrange polynomial foreground basis functions defined upon a boundary-conforming foreground mesh. This work extends the applications of interpolation-based methods to implement RKPM within open-source finite element software. Interpolation-based RKPM is applied to several PDEs, and error convergence rates are equivalent to classic RKPM integrated using high-order Gauss-quadrature schemes. The interpolation-based method is able to exploit the continuity of the RKPM basis to solve higher-order PDEs, demonstrated through the biharmonic problem. The method is extended to multi-material problems through Heaviside enrichment schemes, using local foreground refinement to reduce geometric integration error and achieve high-order accuracy. The computational cost of interpolation-based RKPM is similar to the smoothed gradient nodal integration schemes, offering significant savings over Gauss-quadrature-based meshfree methods while enabling easy implementation within existing finite element software.
In this article, we introduce a semi-orthogonal tribonacci wavelet and develop a semi-orthogonal tribonacci wavelet collocation method, offering an effective numerical method for solving a class of non-linear singular BVPs.
While deep learning has achieved remarkable success in solving partial differential equations (PDEs), it still faces significant challenges, particularly when the PDE solutions have low regularity or singularities. To address these issues, we propose the Weak TransNet (WTN) method, based on a Petrov-Galerkin formulation, for solving elliptic PDEs in this work, though its framework may extend to other classes of equations. Specifically, the neural feature space defined by TransNet (Zhang et al., 2023) is used as the trial space, while the test space is composed of radial basis functions. Since the solution is expressed as a linear combination of trial functions, the coefficients can be determined by minimizing the weak PDE residual via least squares. Thus, this approach could help mitigate the challenges of non-convexity and ill-conditioning that often arise in neural network training. Furthermore, the WTN method is extended to handle problems whose solutions exhibit multiscale features or possess sharp variations. Several numerical experiments are presented to demonstrate the robustness and efficiency of the proposed methods.
This paper proposes an explicit computational method for solving a three-dimensional system of nonlinear elastodynamic sine-Gordon equations subject to appropriate initial and boundary conditions. The time derivative is approximated by interpolation technique whereas the finite element approach is used to approximate the space derivatives. The developed numerical scheme is so-called, high-order explicit computational technique. The new algorithm efficiently treats the time derivative term and provides a suitable time step restriction for stability and convergence. Under this time step limitation, both stability and error estimates of the proposed approach are deeply analyzed using a constructed strong norm. The theoretical studies indicate that the developed approach is temporal second-order convergent and spatially third-order accurate. Some numerical examples are carried out to confirm the theory, to validate the computational efficiency and to demonstrate the practical applicability of the new computational technique.
We present a novel artificial diffusion method to circumvent the instabilities associated with the standard finite element approximation of convection-diffusion equations. Motivated by the micromorphic approach, we introduce an auxiliary variable, which is related to the gradient of the field of interest, and which leads to a coupled problem. Conditions for well-posedness of the resulting formulation are established. We carry out a comprehensive numerical study to compare the proposed methodology against some well-established approaches in one- and two-dimensional settings. The proposed method outperforms established approaches in general in approximating accurately the solutions to pertinent and challenging problems.
We present a general and automated approach for computing model gradients for PDE solvers built on sparse spectral methods, and implement this capability in the widely used open-source Dedalus framework. We apply reverse-mode automatic differentiation to symbolic graph representations of PDEs, efficiently constructing adjoint solvers that retain the speed and memory efficiency of this important class of modern numerical methods. This approach enables users to compute gradients and perform optimization for a wide range of time-dependent and nonlinear systems without writing additional code. The framework supports a broad class of equations, geometries, and boundary conditions, and runs efficiently in parallel using MPI. We demonstrate the flexibility and capabilities of this system using canonical problems from the literature, showing both strong performance and practical utility for a wide variety of inverse problems. By integrating automatic adjoints into a flexible high-level solver, our approach enables researchers to perform gradient-based optimization and sensitivity analyses in spectral simulations with ease and efficiency.
Phase field models have emerged as a powerful and flexible framework for simulating complex interface-driven phenomena across a wide range of scientific and engineering applications. In fracture mechanics, the phase field approach--formulated as a gradient flow of the Griffith fracture energy with Ambrosio-Tortorelli regularization--has gained significant attention for its ability to capture complex crack topologies. In this study, we propose a dynamic fracture phase field model (DF-PFM) based on the elastodynamic wave equation. We further extend this framework by incorporating a unilateral contact condition, yielding a refined model suitable for simulating fault rupture under high pressure. For both models, we rigorously derive energy dissipation identities under mixed boundary conditions, ensuring energy consistency of the formulations. To validate the proposed approach, we conduct numerical experiments using linear implicit time discretization and finite element methods. Our simulations demonstrate that the unilateral contact condition is essential for accurately capturing shear-dominated crack propagation and preventing non-physical interpenetration, especially under high-compression loading scenarios relevant to seismic faulting.
The shallow water equations often assume a constant velocity profile along the vertical axis. However, this assumption does not hold in many practical applications. To better approximate the vertical velocity distribution, models such as the shallow water moment expansion models have been proposed. Nevertheless, under non-slip bottom boundary conditions, both the standard shallow water equation and its moment-enhanced models struggle to accurately capture the vertical velocity profile due to the stiff source terms. In this work, we propose modified shallow water equations and corresponding moment-enhanced models that perform well under both non-slip and slip boundary conditions. The primary difference between the modified and original models lies in the treatment of the source term, which allows our modified moment expansion models to be readily generalized, while maintaining compatibility with our previous analysis on the hyperbolicity of the model. To assess the performance of both the standard and modified moment expansion models, we conduct a comprehensive numerical comparison with the incompressible Navier--Stokes equations -- a comparison that is absent from existing literature.
In this paper, we propose an accelerated version for the Sinkhorn algorithm, which is the reference method for computing the solution to Entropic Optimal Transport. Its main draw-back is the exponential slow-down of convergence as the regularization weakens $\varepsilon \rightarrow 0$. Thanks to spectral insights on the behavior of the Hessian, we propose to mitigate the problem via an original spectral warm-start strategy. This leads to faster convergence compared to the reference method, as also demonstrated in our numerical experiments.
We propose the periodic scaled Korobov kernel (PSKK) method for nonparametric density estimation on $\mathbb{R}^d$. By first wrapping the target density into a periodic version through modulo operation and subsequently applying kernel ridge regression in scaled Korobov spaces, we extend the kernel approach proposed by Kazashi and Nobile (SIAM J. Numer. Anal., 2023) and eliminate its requirement for inherent periodicity of the density function. This key modification enables effective estimation of densities defined on unbounded domains. We establish rigorous mean integrated squared error (MISE) bounds, proving that for densities with smoothness of order $\alpha$ and exponential decay, our method achieves the $\mathcal{O}(M^{-1/(1+1/(2\alpha)+\epsilon)})$ MISE convergence rate with an arbitrarily small $\epsilon>0$. While matching the convergence rate of the previous kernel approach, our approach applies to a broader class of non-periodic distributions. Numerical experiments confirm the theoretical results and demonstrate significant improvement over traditional kernel density estimation in large-sample regimes.
Deep neural networks are powerful tools for solving nonlinear problems in science and engineering, but training highly accurate models becomes challenging as problem complexity increases. Non-convex optimization and numerous hyperparameters to tune make performance improvement difficult, and traditional approaches often prioritize minimizing mean squared error (MSE) while overlooking $L^{\infty}$ error, which is the critical focus in many applications. To address these challenges, we present a progressive framework for training and tuning high-precision neural networks (HiPreNets). Our approach refines a previously explored staged training technique for neural networks that improves an existing fully connected neural network by sequentially learning its prediction residuals using additional networks, leading to improved overall accuracy. We discuss how to take advantage of the structure of the residuals to guide the choice of loss function, number of parameters to use, and ways to introduce adaptive data sampling techniques. We validate our framework's effectiveness through several benchmark problems.
Algorithms for jointly obtaining projection estimates of the density and distribution function of a random variable using the Legendre polynomials are proposed. For these algorithms, a problem of the conditional optimization is solved. Such an optimization allows one increasing the approximation accuracy with a minimum computational costs. The proposed algorithms are tested on examples with different degree of smoothness of the density.