Stochastic gradient methods have gained increasing attention for solving large-scale inverse problems due to their computational efficiency. However, their theoretical justification in the context of ill-posed problems remains underdeveloped, particularly regarding convergence rate analysis, where existing results typically yield only suboptimal rates. In this paper, we address this gap by establishing order-optimal convergence rates for a stochastic gradient method applied to linear ill-posed problems in Hilbert spaces. Under Hölder-type source conditions with smoothness parameter ν∈ (0, 1/2] , we derive convergence rates both in expectation and almost surely, accommodating a broad class of step-size sequences, including constant and polynomially decaying ones. Our analysis is based on a delicate Lyapunov-type argument and an application of the Robbins–Siegmund theorem. As a byproduct, we also establish new convergence results that do not rely on any source conditions.
Hierarchical Model (HiMod) reduction is an effective Reduced Order Modelling technique for problems defined on elongated, pipe-like domains. It is particularly suitable when a dominant dynamics is aligned with the longitudinal direction, while transverse effects are locally significant but spatially limited. When applied to two-field problems such as the Stokes equations, a main challenge is to ensure the stability of the reduced formulation, particularly the inf-sup condition for pressure discretization. In this work, we provide a rigorous analysis showing that the inf-sup condition holds whenever the number of velocity modes is at least equal to the number of pressure modes, thereby extending previous heuristic approaches. The proof exploits the separation of variables in HiMod and is valid for pipe-like domains under some geometric assumptions. Numerical assessment confirms the theoretical findings, providing a solid foundation for stable and efficient HiMod reduction in incompressible flow problems.
We establish fully-discrete a priori and semi-discrete in time a posteriori error estimates for a discontinuous-continuous Galerkin discretization of the wave equation in second order formulation; the resulting method is a Petrov–Galerkin scheme based on piecewise polynomial test functions and continuous piecewise polynomial trial functions in time, respectively. Crucial tools in the a priori analysis for the fully-discrete formulation are the design of suitable projection and interpolation operators extending those used in the parabolic setting, and stability estimates based on a nonstandard choice of the test function; a priori estimates are shown, which are measured in L^∞ -type norms in time. For the semi-discrete in time formulation, we exhibit reliable a posteriori error estimates for the error measured in the L^∞ (L^2) norm with fully explicit constants; to this aim, we design a reconstruction operator into 𝒞^1 piecewise polynomials over the time grid with optimal approximation properties in terms of the polynomial degree distribution and the time steps. Numerical examples illustrate the theoretical findings.
The classical arguments employed when obtaining error estimates of Finite Element (FE) discretisations of elliptic problems lead to more restrictive assumptions on the regularity of the exact solution when applied to non-conforming methods. The so-called minimal regularity estimates available in the literature relax some of these assumptions, but are not truly of minimal regularity , since a data oscillation term appears in the error estimate. Employing an approach based on a smoothing operator, we derive for the first time error estimates for Discontinuous Galerkin (DG) type discretisations of non-linear problems with $$(p,\delta )$$ ( p , δ ) -structure that only assume the natural $$W^{1,p}$$ W 1 , p -regularity of the exact solution, and which do not contain any oscillation terms.
In this article, we study the sparse grid discretization for the numerical solution of the algebraic Riccati equation (ARE). This approach is of particular interest for the solution of large scale AREs. Such AREs arise, for example, from the discretization of operator Riccati equations associated with the linear quadratic control of systems evolving in a Hilbert space H. Following [5, 47], we formulate the ARE as a nonlinear operator equation on the space of Hilbert–Schmidt operators and derive the matrix equation for the sparse grid discretization. If we use N degrees of freedom to discretize the space H, the sparse grid approximation of the ARE has a memory requirement of order 𝒪(N log N) . We further propose an algorithm that evaluates the sparse grid approximation of the ARE with 𝒪(N^3/2) operations. This considerably reduces the cost of solving the ARE compared to the 𝒪(N^2) memory requirement and 𝒪(N^3) complexity of the regular tensor product discretization. Numerical results are presented to validate the approach.
This paper studies the numerical solution of the semiclassical nonlinear Schrödinger equation on the d-dimensional torus 𝕋^d , with highly oscillatory initial data depending on a small parameter ε∈ (0,1] . We first show that a WKB-type approximation attains an 𝒪(ε ) error in the L^2 norm for H^2 initial data theoretically, although its accuracy deteriorates as ε increases. To address this limitation, we propose a numerical scheme that (i) applies a Galilean transform to remove the oscillations in the initial data, (ii) establishes sharp space–time estimates for the transformed equation, and (iii) employs a new low-regularity integrator to achieve second-order accuracy under the minimal H^2 regularity, which is weaker than the regularity assumptions in the literature. Furthermore, our analysis shows that the CFL-type conditions linking h, τ , and ε —typically imposed in the semiclassical regime in the literature—are not required in our scheme to obtain second-order convergence with respect to τ and h, uniformly with respect to ε , under the weaker regularity condition. Numerical experiments support the theoretical results and demonstrate the robustness of the method across a wide range of ε .
We provide the error analysis for one method of orthogonalization of matrix column blocks in floating point arithmetic, which is a crucial step in the one-sided block Jacobi algorithm for computing the singular value decomposition of a general matrix. The orthogonalization is based on computing the Gram matrix and its Cholesky decomposition to obtain the R-factor. Then, the one-sided element-wise Jacobi algorithm is applied to compute the singular value decomposition of the R-factor via Givens rotations, and, finally, the column block is updated by accumulated Givens rotations. We provide the upper bounds for roundoff errors in each of above mentioned steps, and our analysis is based on a row and column scaling of matrices arising during computation both in exact and floating point arithmetic. Our main result is the upper bound for the orthogonality error of computed left singular vectors of a given matrix column block, the form of which is discussed in detail. Numerical experiments illustrate the developed theory.
The finite element approximation of surface evolution under an external velocity field is studied. An artificial tangential motion is designed by using harmonic map heat flow from the initial surface onto the evolving surface. This makes the evolving surface have minimal deformation (up to certain relaxation) from the initial surface and therefore improves the mesh quality upon discretization. By exploiting and utilizing an intrinsic cancellation structure in this formulation and the role played by the relaxation term, convergence of the proposed method in approximating surface evolution in the three-dimensional space is proved for finite elements of degree k≥ 4 . One advantage of the proposed method is that it allows us to prove convergence of numerical approximations by using the normal vector of the computed surface in the numerical scheme, instead of evolution equations of normal vector (as in the literature). Another advantage of the proposed method is that it leads to better mesh quality in some typical examples, and therefore prevents mesh distortion and breakdown of computation. Numerical examples are presented to illustrate the convergence of the proposed method and its advantages in improving the mesh quality of the computed surfaces.
The problem of finding a solution to the linear system Ax = b with certain minimization properties arises in numerous scientific and engineering areas. In the era of big data, the stochastic optimization algorithms become increasingly significant due to their scalability for problems of unprecedented size. This paper focuses on the problem of minimizing a strongly convex function subject to linear constraints. We consider the dual formulation of this problem and adopt the stochastic coordinate descent to solve it. The proposed algorithmic framework, called adaptive stochastic dual coordinate descent, utilizes sampling matrices sampled from user-defined distributions to extract gradient information. Moreover, it employs Polyak’s heavy ball momentum acceleration with adaptive parameters learned through iterations, overcoming the limitation of the heavy ball momentum method that it requires prior knowledge of certain parameters, such as the singular values of a matrix. With these extensions, the framework is able to recover many well-known methods in the context, including the randomized sparse Kaczmarz method, the randomized regularized Kaczmarz method, the linearized Bregman iteration, and a variant of the conjugate gradient (CG) method. Additionally, we introduce an equivalent formulation that, in certain cases, substantially reduces the need for full-dimensional vector operations introduced by the momentum term. We prove that, with strongly admissible objective function, the proposed method converges linearly in expectation. Numerical experiments are provided to confirm our results.
We consider a systematic numerical approximation of a viscoelastic phase separation model that describes the demixing of a polymer-solvent mixture. An unconditionally stable discretisation method is proposed based on a finite element approximation in space and a variational time discretization strategy. The proposed method preserves the energy-dissipation structure of the underlying system exactly and allows us to establish a fully discrete nonlinear stability estimate in natural norms based on the concept of relative energy. These estimates are used to derive order optimal error estimates for the method under minimal smoothness assumptions on the problem data, despite the presence of various strong nonlinearities in the equations. The theoretical results and main properties of the method are illustrated by numerical simulations, which also demonstrate the capability to reproduce the relevant physical effects observed in experiments.
We establish rigorous a posteriori error bounds for a space-time finite element method of arbitrary order discretising linear wave problems in second order formulation. The method combines standard finite elements in space and continuous piecewise polynomials in time with an upwind discontinuous Galerkin-type approximation for the second temporal derivative. The proposed scheme accepts dynamic mesh modification, as required by space-time adaptive algorithms, resulting in a discontinuous temporal discretisation when mesh changes occur. We prove a posteriori error bounds in the L^∞ (L^2) norm, using carefully designed temporal and spatial reconstructions; explicit control on the constants (including the spatial and temporal orders of the method) in those error bounds is shown. The convergence behaviour of the dominant part of the error estimator is verified numerically, also taking into account the effect of the mesh change. A space-time adaptive algorithm is proposed and tested numerically.
Kernel interpolation, especially in the context of Gaussian process emulation, is a widely used technique in surrogate modelling, where the goal is to cheaply approximate an input–output map using a limited number of function evaluations. However, in high-dimensional settings, such methods typically suffer from the curse of dimensionality; the number of required evaluations to achieve a fixed approximation error grows exponentially with the input dimension. To overcome this, a common technique used in high-dimensional approximation methods, such as quasi-Monte Carlo and sparse grids, is to exploit functional anisotropy: the idea that some input dimensions are more ‘sensitive’ than others. In doing so, such methods can significantly reduce the dimension dependence in the error. In this work, we propose a generalisation of sparse grid methods that incorporates a form of anisotropy encoded by the lengthscale parameter in Matérn kernels. We derive error bounds and perform numerical experiments that show that our approach enables effective emulation over arbitrarily high dimensions for functions exhibiting sufficient anisotropy.
We study the discretization of (almost-)Dirac structures using the notion of retraction and discretization maps on manifolds. Additionally, we apply the proposed discretization techniques to obtain numerical integrators for port-Hamiltonian systems and we discuss how to merge the discretization procedure and the constraint algorithm associated to systems of implicit differential equations. After fixing a time step, discretization maps are used to discretize the configuration manifold and we obtain different geometric integrators for the given implicit differential equations. As a result we propose a new method that exactly preserves the integrable part of those equations.
Many-body interactions arise naturally in the perturbative treatment of classical and quantum many-body systems and play a crucial role in the description of condensed matter systems. In the case of three-body interactions, the Axilrod–Teller–Muto (ATM) potential is highly relevant for the quantitative prediction of material properties. This work solves the long-standing issue of the numerical computation of the resulting energies in d-dimensional lattice systems. We present an efficiently computable representation of many-body lattice sums in terms of singular integrals over products of Epstein zeta functions. For three-body interactions in three dimensions, this approach reduces the runtime for computing the ATM lattice sum from weeks to minutes. Our approach further extends to a broad class of n-body lattice sums. We demonstrate that the computational cost of our method only increases linearly with n, evading the exponential increase in complexity of direct summation. We discuss techniques for numerically computing the arising singular integrals and compare the accuracy of our results against computable special cases and against direct summation in low dimensions, achieving full precision for exponents greater than the system dimension. Finally, we apply our method to study the stability of a three-dimensional lattice system with Lennard–Jones two-body interactions under the inclusion of an ATM three-body term at finite pressure, finding a transition from the face-centered-cubic to the body-centered-cubic lattice structure with increasing ATM coupling strength. This work establishes both the numerical and analytical foundation for an ongoing investigation into the influence of many-body interactions on the stability of matter.
We study a variant of the Strang splitting for the time integration of the semilinear wave equation under the finite-energy condition on the torus 𝕋^3 . In the case of a cubic nonlinearity, we show almost second-order convergence in L^2 and almost first-order convergence in H^1 . If the nonlinearity has a quartic form instead, we show analogous convergence results, where the order is reduced by 1/2 in both cases. To our knowledge these are the best convergence results available for the 3D cubic and quartic wave equations under the finite-energy condition. Our approach relies on continuous- and discrete-time Strichartz estimates. We also make use of the integration and summation by parts formulas to exploit cancellations in the error terms. Moreover, error bounds for a full discretization using the Fourier pseudo-spectral method in space are given. Finally, we discuss a numerical example indicating the sharpness of our theoretical results.
This paper investigates quasi-Monte Carlo (QMC) integration of Lebesgue integrable functions with respect to a density function over ℝ^s . We extend the construction-free median QMC rule proposed by Goda and L’Ecuyer (SIAM J Sci Comput, 2022) to the weighted unanchored Sobolev space of functions defined over ℝ^s introduced by Nichols and Kuo (J Complexity 2014). By taking the median of k = 𝒪(log N) independent randomized QMC estimators, we prove that for any ϵ∈ (0,r-1/2] , our method achieves a mean absolute error bound of 𝒪(N^-r+ϵ) , where N is the number of points and r>1/2 is a parameter determined by the function space. This rate matches the rate of randomly shifted lattice rules obtained via a component-by-component (CBC) construction, while our approach requires no specific CBC constructions or prior knowledge of the space’s weight structure. Numerical experiments demonstrate that our method attains an accuracy comparable to the CBC construction based method, and outperforms the Monte Carlo method.
This paper presents a novel approach to rigorously solving initial value problems for semilinear parabolic partial differential equations (PDEs) using fully spectral Fourier–Chebyshev expansions. By reformulating the PDE as a system of nonlinear ordinary differential equations and leveraging Chebyshev series in time, we reduce the problem to a zero-finding task for Fourier–Chebyshev coefficients. A key theoretical contribution is the derivation of an explicit decay estimate for the inverse of the linear part of the PDE, enabling larger time steps. This allows the construction of an approximate inverse for the Fréchet derivative and the application of a Newton–Kantorovich theorem to establish solution existence within explicit error bounds. Building on prior work, our method is extended to more complex partial differential equations, including the 2D Navier–Stokes equations, for which we establish global existence of the solution of the IVP for a given nontrivial initial condition.
This paper studies the solution of nonsymmetric linear systems by preconditioned Krylov methods based on the normal equations, LSQR in particular. On some examples, preconditioned LSQR is seen to produce errors many orders of magnitude larger than classical direct methods; this paper demonstrates that the attainable accuracy of preconditioned LSQR can be greatly improved by applying iterative refinement or restarting when the accuracy stalls. This observation is supported by rigorous backward error analysis. This paper also provides a discussion of the relative merits of GMRES and LSQR for solving nonsymmetric linear systems, demonstrates stability for left-preconditioned LSQR without iterative refinement, and shows that iterative refinement can also improve the accuracy of preconditioned conjugate gradient.
We present a conforming setting for a mixed formulation of linear elasticity with symmetric stress that has normal-normal continuous components across faces of tetrahedral meshes. We provide a stress element for this formulation with 30 degrees of freedom that correspond to standard boundary conditions. The resulting scheme converges quasi-optimally and is locking free. Numerical experiments illustrate the performance.
In this work, we establish that discontinuous Galerkin methods are capable of producing reliable approximations for a broad class of nonlinear variational problems. In particular, we demonstrate that these schemes provide essential flexibility by removing inter-element continuity while also guaranteeing convergent approximations in the quasiconvex case. Notably, quasiconvexity is the weakest form of convexity pertinent to elasticity. Furthermore, we show that in the non-convex case discrete minimisers converge to minimisers of the relaxed problem. In this case, the minimisation problem corresponds to the energy defined by the quasiconvex envelope of the original energy. Our approach covers all discontinuous Galerkin formulations known to converge for convex energies. This work addresses an open challenge in the vectorial calculus of variations: developing and rigorously justifying numerical schemes capable of reliably approximating nonlinear energy minimization problems with potentially singular solutions, which are frequently encountered in materials science.
