Papers

We develop a discrete design theory for cofactor reference forms that preserve the complete finite-element elasticity tangent. A fixed continuous potential generates the coefficient, while the variation space and quadrature determine compatibility. We classify the exact local potential spaces for all-order tensor elements, anisotropic spaces, cubic serendipity and quadratic tetrahedra. At a stress-free Mooney-Rivlin state with a common positive material-coefficient sum, a convex loss expresses the attainable quadrature-point material margin. The optimum is attained; under stated space and sampling assumptions, its ideal value is attained precisely when the physical identity potential is compatible. A conforming fixed-candidate example proves a strict margin loss, linear in curvature amplitude and uniform in cell size, under displacement enrichment. Explicit positive candidates accompany this restriction. A Gram identity gives the sharp constant-coefficient defect at a prescribed margin in the stated norm. Compatible and reference-defect-corrected forms retain the original tangent through a deformation-independent external-boundary Hessian; equal potential traces give equal exact assembled gauges. Quadratic tensor geometry admits a stress-free physical coercivity bound uniform in mesh size and displacement order under explicit material, shape, quadrature and boundary assumptions. Certified curved tetrahedral shear families with zero recovered cell pressure establish finite-deformation regimes, while a loaded pressure-curvature counterexample identifies their limitation. The theory and observations distinguish local exactness, material certificates, complete-reference coercivity and unchanged equations.

Volume quadrature can change nonlinear finite-element constraints while preserving their reference-state derivatives. We connect an explicit determinant defect to feasible-set geometry and singular mechanical response. For affine tensor elements of coordinate degree $p\ge3$ with $n\ge p+1$ Gauss points per coordinate, determinant volume is exact precisely when $2n\ge3p$. Below that threshold we construct a boundary-fixed cubic defect at every order. The same directions yield a full-space cube-root residual–distance bound under explicit cell-support, coefficient and physical-norm assumptions, with mesh-uniform upper constants at fixed order. With all cell-pressure equations retained, the volume Jacobian gains rank at nearby feasible states despite agreement through second derivatives at rest; the cube-root exponent is sharp on each fixed mesh. A general localized-minimum theorem shows that the first reduced compatibility term contributes its weighted square to the leading energy in a joint small-load, large-bulk limit. Cubic and quadratic defects therefore produce sextic and quartic terms. The full cubic-element interior space has an exact normal form and sharp local error exponents. Curved quadratic tetrahedra supply the second-order contrast, a sparse rational witness and an exact four-Jacobian volume formula. Finite-strain tensor calculations illustrate normalized response separation, with explicit stationary-point, extreme-bulk and pressure-recovery qualifications. The constructive correction preserves exact cell volumes, so quadrature feasibility remains distinct from physical volume preservation.

In this work, we recast two widely used acceleration methods for fixed-point iterations, Anderson acceleration (AA) and the nonlinear generalized minimal residual method (NGMRES), as two-grid methods. By explicitly deriving the error propagation matrices for AA and NGMRES on linear problems, we show that both methods, together with several of their existing variants, can be reformulated as two-grid methods with a single pre- or postsmoothing step and a coarse-grid correction given by a projection with respect to a suitable inner product. This reformulation not only unifies AA, NGMRES, and their variants within a single algorithmic framework, but also enables the design of new variants by independently adjusting the components of the two-grid method, such as the coarse space and the smoothing steps. In particular, we propose a new variant that enlarges the coarse space by incorporating the most recently updated residual, an improvement that is naturally revealed by the two-grid formulation but is hidden in the standard formulations of the AA and NGMRES algorithms. Numerical experiments on both an SPD Poisson problem and a family of non-SPD convection-diffusion problems show that this new variant often outperforms both AA and NGMRES, as well as their variants.

We study loss-based filtering for finite-sum optimization with a subset of corrupted component functions whose gradients may be highly unreliable. Motivated by minimum-loss-based SGD (min-$k$-loss) and quantile-based methods for corrupted linear systems, we propose and analyze a general loss-filtering framework – Quantile-\(k\)-Loss SGD (Q\(k\)L-SGD) – that samples \(k\) component losses at each iteration and updates using an index chosen uniformly from the lower empirical \(q\)-quantile. We prove linear convergence of this family of methods under standard convexity assumptions, requiring the sample size to scale with the number of corruptions and a subset strong-convexity threshold. For the cases when large enough sampling is impossible or undesirable, we give a complementary small-sample probabilistic analysis that covers any sample size $k$ and the convergence behavior depends on the probability of selecting an outlier and on the curvature of the selected good step. Experiments on polynomial regression, regularized logistic regression, and regularized hinge loss show that intermediate quantiles often outperform both standard SGD and min-\(k\)-loss SGD. In particular, min-\(k\) often stalls by repeatedly selecting nearly solved components, while intermediate quantiles retain robustness and produce more informative updates.

We develop a unified theory of vectorial quasi-interpolation that exactly preserves intrinsic physical structures characterized by systems of linear constant-coefficient differential operators, including divergence-free and curl-free constraints as special cases. The key ingredient is a family of matrix-valued kernels constructed from an orthogonal projector in the frequency domain, whose columns analytically satisfy the prescribed constraints. The resulting quasi-interpolant is a weighted combination of kernel translates with sampled function values as coefficients, and therefore requires no constrained optimization. Simultaneous error estimates for the approximand and its derivatives are established within a bias-variance framework. For numerical vector-field decomposition, we introduce a generalized Helmholtz–Hodge decomposition that splits a smooth vector field into two components with distinct physical structures, and develop corresponding structure-preserving quasi-interpolation schemes and error estimates. Explicit kernels are derived for several representative constraints, including steady-state acoustic constraints, combined divergence-curl constraints, coupled divergence-free constraints, and second-order Saint-Venant compatibility conditions. Numerical experiments confirm the theoretical error estimates and exact structure preservation, demonstrate accurate identification of sources and sinks, and show three-dimensional reconstruction of generalized Helmholtz–Hodge components for general smooth vector fields without prior constraint assumptions.

We present in this article the discretization of small and finite strain mechanics using variationally consistent, stabilized, nodal strain finite element formulations. We prove that these methods derive from a mixed variational principle that uses different meshes for the primal and dual variables. We demonstrate that, for both mechanical problems, nodal strain finite elements can be written as pure primal formulations with assumed strain operators. The stabilization terms are shown to be necessary and a new, general class of stabilizing functions is proposed for the finite strain regime. The ensuing formulations can be interpreted as pure particle methods because all kinematic and material information is stored at the nodes. This point of view explains the favorable properties of nodal strain methods when employed in problems that require frequent remeshing and involve inelastic materials. In these cases, the diffusion of internal variables is kept to a minimum. The numerical examples presented validate the claims.

In this paper we consider the linear elasticity eigenvalue problem with different boundary conditions and examine two a posteriori error estimators, that compute solutions to local residual problems using a partition of unity and localized residuals. These a posteriori error estimators are known for the source problem, and in this paper we investigate their reliability and efficiency for the given eigenvalue problem with Dirichlet, Neumann and gliding boundary conditions. In numerical experiments we test the performance of the estimators on different domains, with different boundary conditions, and with different polynomial degrees. The numerical experiments show efficiency indices closer to one compared to the standard residual-based a posteriori error estimator.

Multistable systems admit distinct stable equilibria for the same loading and boundary conditions. Shape optimization of multistable structures offers a promising route to engineer the mechanical response of advanced materials and metamaterials. In this work, we develop a formulation and algorithm for controlling the energy jump associated with snap-through behaviour in hyperelastic metamaterials by optimizing the shape of the domain. The key challenge is that classical shape optimization theory assumes a single-valued domain-to-solution map, which does not hold in the multistable setting. We address this by extending the theory to settings with non-unique solutions, treating the scaling of the domain deformation as a continuation parameter and applying the implicit function theorem to give sufficient conditions for the continued existence of multiple solution branches as the shape varies. This theoretical foundation underpins a practical algorithm that targets a prescribed energy jump between stable states, enabling the systematic design of structures whose snap-through response can be tuned on demand.

Higher-order regularization has shown improved results over first-order methods such as total variation. In this paper, convergence of a piecewise linear finite element discretization is shown for the second-order total generalized variation functional. Additionally, this convergence result is extended to a spatiotemporal setting after the spatiotemporal function space is comprehensively introduced. In both cases, discrete minimizers converge to the continuous minimizer with rate $h^{1/4}$, matching the best known rate for finite element discretizations of total variation in the same general setting. Numerical experiments confirm the convergence rates and indicate faster convergence in practice. Comparisons with related spatial image reconstruction algorithms show comparable reconstruction results. The spatiotemporal case is discussed using the example of the inverse problem in electrocardiographic imaging, where the proposed functional improves the state of the art.

This paper presents an optimal interior penalty discontinuous Galerkin (IPDG) scheme for the planar biharmonic equation using piecewise polynomials of degree $k=3$ or $4$. In standard IPDG methods, large penalty parameters force the discrete solution into an overconstrained space, severely degrading accuracy—phenomenon known as numerical locking. To overcome this, our method enforces only vertex continuity and projects the jumps of the function and its normal derivative onto $\mathcal{P}^{k-3}$ and $\mathcal{P}^{k-2}$, respectively. Consequently, as the penalty parameters tend to infinity, the discrete solution is forced to lie in a constrained subspace $V_{h,\infty}^k$, which we identify as the optimal nonconforming finite element space $B_h^k$. This intrinsic connection fundamentally eliminates numerical locking. We prove optimal error estimates of $\mathcal{O}(h^{k-1})$ in a mesh-dependent energy norm and $\mathcal{O}(h^{k+1})$ in the $L^2$ norm. Numerical experiments on both convex and L-shaped domains confirm that the proposed scheme is robust and entirely locking-free, even for extremely large penalty parameters.

Multi-electrode Arrays (MEAs) enable tissue-level electrophysiological studies of human induced pluripotent stem cell-derived cardiomyocytes (hiPSC-CMs) by recording extracellular field potentials. Recent computational models have improved MEA simulations by coupling detailed electrode descriptions with the Bidomain framework, a parabolic-elliptic system of nonlinear PDEs coupled with a stiff ionic model. The standard numerical strategy relies on operator splitting techniques that decouple the PDE and ODE components. In most implementations, the ionic subsystem is treated explicitly while the diffusive operator is handled implicitly, resulting in a first-order implicit-explicit (IMEX) time discretization. Although computationally convenient, this approach limits temporal accuracy and may reduce efficiency in large-scale simulations. In this work, we investigate higher-order IMEX Runge-Kutta schemes within a Strang-based operator splitting, specifically tailored for the MEA model of hiPSC-CMs monolayers. First- and second-order schemes are compared in terms of computational cost and global error against a high-fidelity reference solution obtained with a very small time step. For a fixed computational cost, the second-order schemes achieve errors that are 2-3 orders of magnitude smaller than those of the first-order method. These results demonstrate that higher-order IMEX integration significantly improves the accuracy-to-cost ratio of MEA simulations, providing a practical and reliable approach for large-scale electrophysiological studies.

We present a biparametric preconditioning technique for large, sparse, non-symmetric, and non-singular double saddle point problems. The preconditioner is induced using a stationary iteration method. It is also based on a two-parameterized block LU factorization of the coefficient matrix. To begin with, some convergence results of the iteration method are derived. Moreover, the spectral features exhibited by the preconditioned matrix are given. Finally, the presented preconditioner is conjuncted with the GMRES method. Numerical results demonstrate the satisfactory performance of the preconditioner.

We analyse fully discrete Streamline Diffusion finite element methods, also known as SUPG methods, for time-dependent first-order systems with skew-symmetric spatial operators. To the best of our knowledge, this is the first stability and a priori error analysis of strongly consistent Streamline Diffusion finite elements combined with linear multistep schemes of order greater than two for multidimensional first-order hyperbolic systems. Local-in-time $L^2$-residual minimisation singles out $\delta=b_0\tau$, where $\delta$ is the stabilisation parameter, $\tau$ the time step, and $b_0$ the current-time coefficient in the method average. For this choice, each implicit linear system is a symmetric positive definite normal equation in a graph norm. A direct dual-residual estimate yields errors of order $O(h^{k+1/2}+\tau^\nu)$ when $\delta$ is proportional to $h$ and $\tau=O(h)$, where $h$ is the mesh size, $k$ the polynomial degree, and $\nu$ the temporal order. Stability is established for the $\theta$-method and the Adams-Moulton methods AM3-AM5. The general-degree AM3 and AM4 estimates use strengthened CFL conditions; for continuous piecewise affine elements, an elementwise cancellation recovers the standard hyperbolic regime. AM5 and the explicit Adams-Bashforth methods AB3 and AB4 are stable under a standard hyperbolic CFL condition, though the explicit schemes do not have normal-equation structure. Order-matched acoustic tests recover orders two through five and material-residual rates one half-order lower. In a d'Alembert test with discontinuous initial data, normal-equation SUPG recovers second-order local $\mathbb P_1$ convergence away from the wave fronts, whereas the unstabilised error is close to half order. A two-dimensional compact-wave test shows improved localisation with modest dissipation.

This article establishes sharp strong error estimates for the fully discrete approximation of semilinear parabolic stochastic partial differential equations (SPDEs) driven by multiplicative noise and subject to nonsmooth initial data. The spatial discretization is based on a standard finite element method, coupled with the linearly implicit Euler scheme in time. By mapping the diffusion operator into negative fractional spaces, our framework accommodates both trace-class and space-time white noise. For nonsmooth initial data, by decoupling the noise regularity parameter $\beta \in (0,2)$ from the initial data regularity parameter $\mu \in (0,2]$, we derive sharp regularity estimates that isolate the exact loss of initial regularity into an integrable temporal singularity. Furthermore, we establish sharp strong convergence rates of $O(h^{\beta-\varepsilon} + k^{\frac{1}{2}\min\{\beta-\varepsilon, 1\}})$ for $\varepsilon>0$ away from $t = 0$. Finally, we consider physically relevant stochastic models, such as the modified Langmuir fractional surface coverage model and the parabolic Anderson model, in our numerical experiments to confirm the theoretical convergence rates.

In this paper, we study the structural relationship between immersed finite element (IFE) spaces for scalar elliptic and Stokes interface problems on unfitted meshes. Although IFE methods have been widely developed for individual partial differential equation models, the algebraic connection between scalar and vector IFE spaces has not been systematically explored. We establish a precise unisolvence relationship between the immersed Crouzeix-Raviart (CR) element for elliptic interface problems and the immersed CR-$P_0$ element for Stokes interface problems in both gradient and stress formulations. By analyzing the block structure of the local IFE matrices, we show that the determinant of the Stokes IFE matrix factorizes in terms of the determinant of the corresponding elliptic IFE matrix. This factorization reveals an intrinsic algebraic link between the two classes of immersed elements and explains the unisolvence of the Stokes IFE spaces through the scalar elliptic case. The same result extends naturally to three dimensions and significantly simplifies the corresponding analysis. Numerical experiments confirm optimal convergence rates for both velocity and pressure in three-dimensional Stokes interface problems.

We study an online sketched Newton method that approximates the Newton direction at each step via a state-of-the-art sketching solver, called the generalized accelerated sketch-and-project solver (GAS), thereby mitigating the computational bottleneck of classical second-order methods. The GAS solver improves upon vanilla, unaccelerated sketch-and-project solvers by achieving accelerated convergence through Nesterov momentum updates, and accommodates a flexible projection metric whose proper choice further reduces computational cost. Building on this design, we establish asymptotic normality of the averaged sketched Newton iterates and characterize their limiting covariance matrix. The resulting covariance recovers that of the unaccelerated sketched Newton method under a specific choice of acceleration parameters, converges more rapidly (in the number of sketching steps) to the minimax-optimal covariance in general, and is smaller than that of the last iterate produced by the accelerated method. Finally, we strengthen these results by establishing a functional central limit theorem for the Newton iterates, which allows us to bypass explicit covariance estimation and develop an online inference procedure based on random scaling. Specifically, we construct a pivotal test statistic by appropriately rescaling the averaged iterates, so that its limiting distribution is free of any unknown parameters, enabling asymptotically valid online inference. Numerical experiments demonstrate superior performance of the proposed inference procedure.

Microwave ablation (MWA) is a percutaneous ablation that induces thermal necrosis by generating localized high-temperature regions within tumor tissue. However, when a tumor is located close to a large vessel such as the aorta, blood flow continuously removes heat from the ablation zone. This phenomenon causes the vascular heat-sink effect and thereby increases the risk of incomplete ablation. This study constructs an ANSYS-based numerical framework to investigate the heat-sink effect during MWA of tumors located close to the aorta. The computational model consists of the tumor, surrounding tissue, vessel wall, blood flow domain, and microwave antenna. Microwave-induced heating was represented by an equivalent heat source mapped from the High-Frequency Structure Simulator (HFSS)-derived volume loss density, rather than by direct electromagnetic-thermal coupling within the thermal-fluid solver. Two parametric studies were performed to quantify the heat-sink effect. First, the tumor-aorta distance was varied at D = 0 mm, 4 mm, and 8.5 mm to assess how vessel proximity affects temperature elevation and ablation zone formation. Second, at each tumor-aorta distance, three blood-flow conditions were compared: no blood, stagnant blood V = 0 m/s, and normal aortic blood flow V=0.65 m/s. This comparison enabled the heat-sink effect of blood flow to be independently evaluated. The results highlight the clinical need for simulation-based treatment planning in patients with tumors located close to major vessels.

Stochastic variance reduced gradient (SVRG) is a variant of stochastic gradient descent and is a promising iterative method for solving large-scale inverse problems. Nevertheless, the development of theoretically grounded a posteriori stopping rules for SVRG remains an open challenge. In this work, we provide a convergence analysis of SVRG equipped with the discrepancy principle, the most well-known a posteriori stopping rule, for solving a class of linear inverse problems in Hilbert spaces. We establish the regularizing property of SVRG, and moreover, under suitable source conditions, we derive convergence rates of SVRG iterates. To the best of our knowledge, these are the first convergence rate results of any stochastic iterative method for inverse problems under the a posteriori stopping rule. The theoretical findings are supported by numerical experiments.

We propose spectral and spectral-element methods for accurately computing many eigenvalues of the Schrödinger operator with Robin boundary condition on simple and complex geometries, respectively. Due to their spectral-type accuracy in approximating those eigenfunctions corresponding to high-index eigenvalues, which are usually highly oscillatory, the proposed approaches have excellent resolution in computing thousands of eigenvalues accurately and efficiently with the resolution property that the number of eigenvalues with reasonable accuracy is proportional to the number of degrees of freedom. Using a standard HPC node with 128 GB of memory, we can obtain numerically more than 5,000 reliable eigenvalues with relative errors below $10^{-8}$ for different complex domains in two dimensions (2D). Based on these computed eigenvalues, we first confirm some theoretical results on the Robin-to-Neumann (RtN) gaps of the Laplacian operator in 2D, which were recently studied by Rudnick et al. [Comm. Math. Phys. 388 (2021)]. Then we systematically study the RtN gaps of the Schrödinger operator and their convergence rates. Based on our extensive numerical results, we formulate a unified conjecture on the cumulative averages of the RtN gaps of the Schrödinger operator.

Motivated by accelerating machine learning (ML), many computer vendors and chip manufacturers are building accelerators for matrix multiplication, which save time and energy by operating in the lower precisions needed for ML. This has in turn motivated many efforts to use these accelerators to provide faster matrix multiplication implementations with the higher precision required by many other linear algebra applications. Motivated by the large design space of algorithms for approximating higher precision, with significant performance/accuracy tradeoffs, we provide a benchmark to "grade the accuracy" of a matrix multiplication implementation (or the BLAS more generally), ranging from an "A" for attaining the classic floating point error bound, to a "C" for attaining a weaker but still useful bound, that is satisfied by Strassen-like algorithms. We also propose "ungameable" tests that vendors or users can run to verify their promised accuracy, and describe how these different grades impact the accuracy of applications like LU, QR, and Cholesky decomposition. Our test code is publicly available at github.com/Reference-LAPACK/grading-the-BLAS for developers and users, and we also plan to publicly release all our test results.

Many scientific applications require modeling how diffusive systems evolve over time, not merely their eventual steady states. While conventional steady-state analysis of partial differential equations (PDEs) on complex geometries is already hindered by costly volumetric meshing, transient analysis further requires sequential time stepping and careful step size selection. Grid-free Monte Carlo solvers such as walk on spheres (WoS) and walk on stars (WoSt) avoid this meshing bottleneck but remain largely limited to steady-state problems. We generalize WoS, for pure Dirichlet problems, and WoSt, for mixed Dirichlet–Neumann problems, to heat equations with initial conditions and time-dependent source and boundary data. We equip each random walk with a finite time budget and sample an exit time at every spatial step. If the exit time exceeds the remaining budget, the walk samples an interior point and evaluates the initial condition; otherwise, it continues with a reduced budget, accumulating source and boundary contributions. Our main technical contribution is a suite of kernel sampling and variance reduction techniques, including a low-bias, tabulation-free exit time sampler and efficient rejection samplers. Unlike grid-based transient solvers, our method directly estimates the solution at any requested time without volumetric meshing or sequential time marching. It also retains the parallel, progressive, and output-sensitive evaluation of WoS and WoSt while eliminating time step selection and temporal discretization bias entirely. Finally, we show how sharing walks enables efficient estimates at multiple target times.

We study the selection and stable realization of local approximation spaces in Trefftz discontinuous Galerkin discretizations of the Helmholtz equation. A scaled Cauchy-trace inner product places plane waves and Fourier-Bessel functions in a common geometry: Fourier-Bessel modes are orthogonal with explicit weights, while the same weights determine the circulant spectrum of an equispaced plane-wave trace Gram matrix. This separates amplitude scaling from genuine trace-rank loss and yields an exact best-approximation identity for mixed plane-wave–Fourier–Bessel spaces. With complex plane-wave angles, the unresolved modal tail is an exponential sequence, so propagating and evanescent components can be identified by the same ESPRIT/variable-projection procedure. We prove exact recovery and a perturbation estimate for the recovered angles and, under the standard PWDG quasi-optimality bound, transfer these perturbations to the DG error. Trace-Riesz orthonormalization is then separated from a graph–Riesz normalization of the assembled operator. Numerical experiments verify the identities, recover sparse ray fields to roundoff, and resolve a propagating-to-evanescent transition without a prescribed critical angle.

It is known that the Sherman–Morrison (SM) formula is not numerically stable. In recent work, we introduced SMIR, an algorithm that incorporates iterative refinement to enhance the SM backward error. In this paper we take a different route: adapting an algorithm of Govaerts, originally designed for general bordered linear systems, we develop a modified Sherman–Morrison (MSM) method whose built-in self-correction makes it surprisingly resilient. Whereas standard SM requires the solution of two $n\times n$ linear systems, MSM requires three; by contrast, SMIR requires $2+k$ solves, where $k$ is the number of IR steps and can be substantially larger than three, when IR converges slowly. We then derive backward and forward error bounds for both SM and MSM. The SM backward error bound established here is stronger than the one proved in [Hashemi & Nakatsukasa 2026]: it accounts for every rounding error and holds with no conditions on the size of the capacitance. From these bounds we extract growth factors that are cheap to compute a posteriori and can be used to certify the backward and forward stability of SM and MSM on a given problem. In our experiments, MSM consistently produces backward stable solutions, and is therefore observed to be forward stable as well. A formal proof of the stability — or instability — of MSM remains an open problem.

For $1<p<2$, the $p$-Laplace equation and its $p$-Stokes generalization are difficult to solve numerically. Newton's method converges rapidly only close to the solution, with iteration counts that grow under mesh refinement and deteriorate as $p\to1$. The more robust Picard iteration converges only linearly. Rather than globalizing or preconditioning Newton's method, we modify the system to which it is applied. We lift the equation by introducing the flux $|\nabla u|^{p-2}\nabla u$ as an auxiliary (or ”lifting”) variable, apply a nonlinear transformation to the resulting constitutive relation, linearize, and eliminate the auxiliary variable by static condensation. Lifting alone leaves the linearization unchanged; it is the preceding transformation that yields the new method. The elimination is algebraic and pointwise at the quadrature points, so the flux variable is never discretized, no inf-sup condition or indefinite system arises, and the cost per iteration is that of a standard Newton step. Together with a pointwise feasibility bound on the lifting variable that keeps the diffusion tensor uniformly positive definite, this yields an iteration that we prove, in finite dimensions and for the $p$-Laplace equation, to converge globally and locally at a quadratic rate; a one-dimensional model problem explains why the lagged flux variable removes the zig-zag behavior of standard Newton for $p$ close to one. Firedrake-based experiments for $p$-Laplace problems in two and three dimensions and for stationary and time-dependent $p$-Stokes flows show iteration counts largely insensitive to $p$ and to mesh refinement, and up to an order of magnitude fewer iterations than standard Newton for $p$ close to one.

In this work we consider numerical models with random input parameters, where the model output is a spatially distributed random field. We carry out a statistical sensitivity analysis of spatial sets which arise for instance when the model output exceeds a critical threshold. We consider two approaches: (i) a kernel-based sensitivity analysis working with the Hilbert–Schmidt Independence Criterion, and (ii) a function-valued sensitivity analysis working with generalized Sobol' indices, where the model output is the indicator function of the spatial set of interest. We develop efficient estimators for the sensitivity indices, especially for finite-element-based numerical models, where we approximate volume integrals in the sensitivity measures by finite element quadrature, thereby avoiding Monte Carlo sampling over the spatial domain. This improves the state-of-the-art for the kernel-based indices and makes these methods practically accessible for expensive numerical models. For the generalized Sobol' indices we utilize the same finite element quadrature for efficient computation. We present numerical experiments for a hydrogen combustion process, modeled by a coupled system of convection-diffusion-reaction equations in a two-dimensional spatial domain, and ask whether the temperature remains below a critical value in selected regions of the combustion domain. Our results show that in this example the kernel-based and the generalized Sobol' indices give qualitatively similar input importance rankings. Moreover, the ranking depends on the chosen observation window and can flip between different windows. The algorithms for implementing the proposed methods are provided in the supplementary materials and in an open source code repository.

In this paper, we develop and analyze a fully discrete local discontinuous Galerkin (LDG) method with IMEX-Euler time discretization for a class of multi-dimensional quasilinear stochastic convection-diffusion-type equations driven by multiplicative $\mathcal Q$-Wiener noise. The leading diffusion matrix may depend on the solution as well as the spatial and temporal variables, while the lower-order drift and noise coefficients may depend on both the solution and its gradient. Under a suitable stochastic parabolicity condition, we establish unconditional high-moment stability estimates for the fully discrete scheme in the quasilinear setting. In the semilinear setting, where the leading diffusion matrix is independent of the solution but may vary in space and time, we further prove optimal high-moment strong error estimates of order $\mathcal O(h^{r+1})$ in space and $\mathcal O(k^{1/2})$ in time. A pathwise error estimate is then derived by combining the high-moment error bound with a discrete Kolmogorov argument. Numerical experiments are presented to illustrate the stability and convergence properties of the proposed method.

We study limited memory steepest descent (LMSD) with exact algebraic rank compression for strictly convex quadratic optimization. At each cycle, the method restricts the gradient history to its column space and applies the reciprocals of all Ritz values of the compressed projection in the next sweep. If the block-start gradient has at most $p$ active distinct eigenvalues, the delayed sweep terminates finitely. For nonterminating trajectories, a determinant and Cauchy–Binet formula for the complete-sweep polynomial yields global convergence without a full-rank history or a run-wise normalized-conditioning bound. Consecutive sweep endpoints define a continuous positively homogeneous map on a closed cone of compatible states. Compactness then gives an $R$-linear endpoint estimate, uniform over compatible initial states, with constants depending only on $H$ and $p$; the estimate extends to inner gradients, iterate errors, and objective gaps. Every fixed positive spectral weight $W=\omega(H)$ is linearly conjugate to the standard method. The weighted methods therefore share its decay factor, while the Euclidean-norm prefactor is increased by at most $\sqrt{\kappa_2(W)}$. This includes harmonic-Ritz LMSD, fixed power weights, and the delayed BB1 and BB2 recurrences.

We give an $R$-linear convergence proof for BB1, BB2, and fixed positive spectral-weight variants on finite-dimensional strictly convex quadratics. The delayed recurrence is written as a first-order system on the closed cone of compatible consecutive-gradient pairs. Its transition is continuous and positively homogeneous of degree one, including at finite termination. A compactness theorem shows that pointwise convergence of every cone orbit implies a uniform finite-step contraction and hence a global $R$-linear estimate. By a realization lemma, every compatible state initiates a BB1 trajectory, so Raydan's theorem yields pointwise stability on the entire cone. Finally, the transformation $g\mapsto W^{1/2}g$ conjugates every fixed positive spectral-weight rule to BB1. Consequently, BB1, BB2, and all such weighted rules have the same homogeneous growth radius.

The Barzilai–Borwein (BB) method is $R$-linearly convergent on strongly convex quadratics, although neither the objective value nor the gradient norm is generally monotone. We construct an explicit current-point Lyapunov function for the quadratic dynamics after the warm-up step, assuming that both endpoints of the initial active spectrum survive this step. If $a<b$ are these endpoints and $P_a,P_b$ are the corresponding spectral projectors, then \[ \overline{f}(x) :=\lVert P_a\nabla f(x)\rVert^{\frac{2b}{a+b}} \lVert P_b\nabla f(x)\rVert^{\frac{2a}{a+b}} \] satisfies, for BB1, BB2, and every fixed positive weighted delayed Rayleigh rule, \[ \overline{f}(x_{k+1}) =\left(\frac{b-a}{b+a}\right)^2 e^{-2D_{a,b}(\tau_k)}\overline{f}(x_k), \qquad D_{a,b}(\tau_k)\ge0, \] where $\tau_k=1/\alpha_k$ is the reciprocal step. The function is obtained by denormalizing the endpoint coboundary certificate in the sharp-rate analysis of Yang and Yuan. We prove that the endpoint exponents are the unique minimax choice among endpoint monomials and illustrate the Lyapunov law on a nonmonotone BB1 trajectory. Using this law, we also give a proof of $R$-linear convergence of the method in finite dimensions.

We propose a learning procedure for system identification in interacting particle systems from single-snapshot observations of collective behaviors, unlike existing approaches that rely on observations of trajectories. This setting leads to a fundamentally ill-posed inverse problem, which we solve by using a regularization strategy based on the empirical distribution of observed configurations, drawn from different, unobserved initial conditions. We test our learning procedure on a variety of representative models with steady-state and quasi-stationary patterns, where collective behaviors encode implicit information about the interaction mechanisms, demonstrating that our approach enables stable and accurate recovery of the underlying interaction laws, leading to faithful reproduction of the collective behavior, and in many cases even of the dynamics leading up to it.