Papers

This article develops optimal order (up to logarithmic factors) weighted cubature rules for Monte Carlo numerical integration in Sobolev reproducing kernel Hilbert spaces over common matrix manifolds. Specifically, we consider matrix models of the Stiefel, Grassmann, flag, and Lagrangian Grassmann manifolds each over either real numbers $\mathbb{R}$, complex numbers $\mathbb{C}$, or quaternions $\mathbb{H}$. The proposed cubature rules use evaluation points which are independent and identically distributed uniform on the manifold. The weights are chosen to be optimal with respect to the reproducing kernel Hilbert space worst-case error, under the constraint that the weights sum to one. Crucially, this constraint ensures optimal weights may be found without having to evaluate the kernel mean embedding (the integral of the kernel), which is often difficult to compute. We consider computable Matérn kernels of Frobenius distances between matrices which reproduce the Sobolev reproducing kernel Hilbert spaces with equivalent norms. Open-source numerical experiments validate our theoretical development, including an application to simultaneous approximation of intrinsic volumes in the Cauchy–Kubota formula. We also present preliminary empirical results beyond independent points including transformations of low-discrepancy quasi-Monte Carlo points from the unit cube and jointly optimized points and weights.

Many dense matrices in scientific computing have hierarchically low-rank off-diagonal structure but are accessible only through routines that apply the matrix and, in general, its transpose or adjoint to vectors, that is, through matrix–vector products (matvecs). We present a fully matrix-free algorithm that compresses such a matrix into Hierarchically Semi-Separable (HSS) form and factorizes it simultaneously with the off-diagonal ranks chosen adaptively from a prescribed tolerance. The method uses $\mathcal{O}(k)$ matvecs with each operator, where $k$ is the largest off-diagonal rank found. Building on fixed-rank block nullification, it samples the matrix in small random batches and enlarges each local basis until a randomized error estimate falls below the tolerance. Unitary nullspace updates give $\mathcal{O}(k^3)$ local dense work, the same asymptotic cost as when the ranks are known in advance. Because compression and elimination proceed together, no separate factorization step is required. Excluding the matvecs, the factorization uses $\mathcal{O}(N)$ arithmetic and storage when the off-diagonal ranks are bounded, and its inverse is applied in $\mathcal{O}(N)$ time. Experiments on Poisson and Helmholtz Schur complements, Matérn covariance matrices, and a complex Lippmann–Schwinger equation demonstrate the factorization as both a direct solver and a preconditioner on problems with up to millions of unknowns.

We establish quantitative weak convergence for the joint terminal law of the solution to a stochastic wave equation and the normalized error of its exponential Euler approximation. For multiplicative trace-class noise and spatial regularity $\alpha\in(1,2]$, the normalized error converges to its limiting distribution at rate $O(m^{-(\alpha-1)})$, where $m$ is the number of time steps. This estimate identifies the leading weak bias and yields a first-order weak-error expansion with remainder $O(m^{-\alpha})$. Moreover, by canceling the leading weak-error term, we construct a Richardson extrapolation scheme with weak convergence order $\alpha$ for $\alpha\in(1,2]$ and show that this rate is sharp. To overcome the lack of parabolic smoothing and control the remainders amplified by normalization, we combine dimension-uniform regularity estimates for the joint Kolmogorov equation with a conditional weak expansion and backward transport of coefficient directions.

We derive explicit formulas for the dual Moore–Penrose generalized inverse (DMPGI) of dual complex matrix sums under some conditions. We first express the DMPGI of a matrix partitioned into column blocks in terms of these blocks and associated projectors, and then derive two formulas for matrix sums. Two examples illustrate the results. We further apply the formulas to rigid-body motion estimation with two observation points subject to a first-order perturbation of their positions that preserves the rank-deficient observation geometry. The DMPGI recovers the observable motion components from the velocity data.

We consider newly designed matrix algebras associated with Cauchy orthogonal transforms. In this work, we are interested in checking the approximation of Toeplitz matrices using algebras of matrices, where each matrix is diagonalized by a fixed Cauchy orthogonal transform, in the spirit of fast algorithm based on preconditioning. We use the Korovkin theory for checking the Toeplitz approximability. Unfortunately and somehow surprisingly, the Korovkin test is negative, in contrast with the successful test for all the Fourier-like and trigonometric algebras, widely used in the relevant literature. More in detail, for the Cauchy algebras, the Frobenius optimal approximation fails to lead to the desired clustering, as also confirmed by ad hoc numerical experiments. It remains open to understand which kind of structured matrices can be approximated by these new algebras.

The latent variable proximal point (LVPP) method solves optimization problems in function spaces with pointwise constraints by a sequence of smooth subproblems whose iterates are strictly feasible. We apply it to linear-quadratic elliptic optimal control problems. Control and observation are distributed but may be restricted to subdomains. We prove that the LVPP subproblems are well posed in appropriate function spaces, that the iteration converges for every sequence of proximity parameters with divergent sum, and that it delivers the Lagrange multiplier of the constraint as a by-product. We interpret the Newton step for the LVPP subproblem as a linear-quadratic optimal control problem with a spatially varying control cost, and we show that Newton's method converges locally quadratically. The Newton systems are solved by the minimal residual method with a block-diagonal preconditioner. We prove that the spectrum of the Newton operator with respect to the inner product defined by the preconditioner is bounded independently of the proximity parameter, the iterate, the mesh and the observation and control domains. The results are combined into a complete algorithm, and tested in numerical experiments.

The sums of the absolute values of higher-order cumulants are classical contrast functions for independent component analysis (ICA) and lead naturally to robust approximate diagonalization of higher-order cumulant tensors using $\ell_1$-based diagonalization over the orthogonal or unitary group. We propose Jacobi-SubG, a Jacobi-type framework that selects index pairs using Riemannian subdifferential information and performs plane rotations with a refinement step that preserves nonzero diagonal contractions. For the general tensor model, we prove global convergence of the entire iterate sequence from any initialization with nonzero diagonal contractions, provided that the subproblem solvers satisfy a uniform sufficient-ascent estimate. For real and complex third- and fourth-order cumulant tensors, which are widely used in ICA, we develop specialized subproblem solvers satisfying these conditions. The resulting algorithms thus enjoy global convergence under the simple initialization condition that all diagonal contractions are nonzero. Numerical experiments on ICA and tensor diagonalization show that the proposed Jacobi-SubG variants achieve better source recovery and higher final objective values in several settings, particularly at the higher tested noise levels.

This work presents a hybrid numerical approach for solving linear systems arising from the discretization of the two-dimensional parametric Helmholtz equation. A convolutional neural network based on the U-Net architecture is trained in an self-supervised manner (i.e., without explicitly generating reference solutions) to approximate the inverse of the discretized Helmholtz operator, using a loss function involving the residual norm of the linear system. The trained network is used as a nonlinear preconditioner within the flexible GMRES (FGMRES) algorithm. Numerical experiments show that while the neural network is not accurate enough to act as a standalone solver, it significantly improves the convergence of FGMRES when employed as a preconditioner. The neural preconditioner demonstrates robust performance and generalization capabilities with respect to variations in the velocity field and the domain size. Comparisons with algebraic preconditioners, based on LU factorizations of approximations of the original system matrices, indicate the superior efficiency of the neural approach under equivalent conditions. We believe that the proposed method is not tied to a specific neural architecture and can be extended to other parametric PDEs.

We prove that the restarted Arnoldi algorithm for real symmetric matrices with restart length $s=2$ has a cyclic limiting behavior. In the terminology established in (arXiv:2209.14579), this is a proof of the ACI($s$) conjecture for $A=A^T\in \mathbb{R}^{n\times n}$ and $s=2$. As explained in (arXiv:2209.14579), this also proves the Forsythe conjecture (Numerische Mathematik, 11 (1968), pp. 57-76) about the limiting behavior of the conjugate gradient method and positive definite $A=A^T\in\mathbb{R}^{n\times n}$ for $s=2$. The proof we give is an extensively rewritten and simplified version of a proof found using Codex with GPT-5.6 Sol. At the end of this paper, we discuss its origin and evolution, as well as related recent developments around the Forsythe conjecture. Our motivation is to present a proof that is instructive and transparent, so that it can enhance the knowledge of interested human readers.

The robustness of two-level overlapping additive Schwarz (OAS) preconditioners for highly heterogeneous elliptic problems hinges on the choice of the coarse space. We propose an algebraic multiscale (AM) coarse space built from the system matrix alone through an algebraic formulation of the multiscale finite element method (MsFEM). The coarse basis functions are energy-minimizing extensions of the solutions of reduced boundary condition (RBC) problems on the subdomain interface, so that coefficient variations along the interface enter the coarse space. By posing the RBC problems on slabs that extend the interface, the construction also applies to general unstructured domain decompositions. Embedding the AM formulation in an energy-minimizing framework extends it to linear elasticity, and we show that existing condition number bounds for OAS with MsFEM coarse spaces carry over to the algebraic setting. For problems with high coefficient inclusions on the coarse nodes, the AM iteration count and condition number are unaffected by increasing the contrast. In 3D scalar diffusion with log-normal coefficients, AM needs about as many iterations as the generalized Dryja–Smith–Widlund (GDSW) coarse space, while solving a coarse problem ten times smaller. For linear elasticity, omitting the coarse basis functions associated with the rotational rigid body modes yields a fully algebraic coarse space that still scales well for homogeneous problems.

We revisit the problem of approximating a linear time-invariant LTI model by a reduced-order model (ROM) in an $H_2$-optimal sense. A widely used method for computing such ROMs is the iterative rational Krylov algorithm (IRKA). However, in practice, the stability of the resulting ROM is not generally guaranteed. In this work, we develop a novel framework that can be used to obtain $H_2$-optimal models of single-input single-output LTI systems while enforcing stability of the ROM. We use Malmquist-Takenaka orthogonal expansions in the corresponding $H_2$-spaces and their connection to hyperbolic geometry to develop an equilibrium function. We show that the poles of an $H_2$-optimal ROM transfer function are defined by the zeros of this equilibrium function. The resulting nonlinear root-finding problem is solved using a Newton-type method formulated on the products of Poincaré disk manifolds, where updates are performed along hyperbolic geodesics ensuring that the stability constraint is satisfied at every iteration. We present numerical experiments on benchmark systems for both discrete- and continuous-time LTI MOR problems, including examples where IRKA produces unstable reduced models while the proposed method preserves stability throughout the iteration.

A ghost-cell immersed boundary method is proposed for the two-dimensional shallow water equations with a variable bed (topography), solved with an explicit Godunov-type finite volume scheme (unsplit second-order MUSCL, HLLC flux, hydrostatic reconstruction) on Cartesian grids. Ghost-cell methods leak mass through the wall. To find a robust and conservative treatment, eight constructions of the ghost values and a cut-cell form of the continuity equation were applied to seven test problems, which include wetting and drying on sloping beds and a laboratory dam break. In the resulting scheme, the ghost values are interpolated bilinearly at the image point; solid and dry nodes are excluded, and the weights of the remaining nodes are renormalized. Mass is conserved by solving the continuity equation in cut-cell form, while the momentum equations are left unchanged. The proposed scheme is about as accurate as the ghost-cell method without this modification and preserves a lake at rest whose shoreline crosses an immersed wall. Because only the continuity equation is modified, it avoids the treatment of the momentum equations on cut cells (wall pressure, cell merging) that a full cut-cell method requires.

This article develops matrix-free Q-less QR algorithms for orthogonalizing tall linear operators. Our main algorithm, CQRRTO, uses randomized preconditioning to make a single Cholesky QR pass numerically safe for any numerically full-rank operator, overcoming the stability issues of Cholesky QR and related methods that operate on the Gram matrix. Compared with deterministic stabilized variants such as Q-less Cholesky QR2 (reorthogonalized Cholesky QR), it remains numerically stable under substantially weaker conditioning assumptions while requiring fewer operator applications. Matrix-free Q-less QR can be useful for improving the numerical robustness of least-squares solves via the normal equations, spectral approximations in a low-dimensional basis, and inversion of operators with low-rank updates via the Woodbury identity. We present a componentwise stability analysis of preconditioned Q-less Cholesky QR that applies to matrix-free operators and is modular with respect to the preconditioner, yielding guarantees for Q-less Cholesky QR, Cholesky QR2, shifted Cholesky QR3, and CQRRTO as special cases. All methods are implemented in the open-source RandLAPACK library.

We use centered differences to construct Chernoff families and derive even asymptotic expansions for their iterated products. Richardson extrapolation is then used to cancel polynomial part and thereby raise the convergence order. As illustrations, we apply the construction to uniformly parabolic, Schrödinger-type second-order operators and to a third-order dispersive operator, obtaining, at extrapolation depth $m$, the rates $ O(n^{-2m}) $ on the Sobolev scales $H^{2m+\delta}$, $H^{4m}$, and $H^{5m}$, respectively, with arbitrary $\delta>0$. In the analytic category, quantitative all-order estimates allow the extrapolation depth to grow with the base resolution and yield the super-fast rates $e^{-cn^2}$ for the two second-order models and $e^{-cn} $ for the third-order model. The constructions are explicit in terms of the coefficients of the underlying differential operators and do not use the target semigroup in their definition. This provides a constructive positive answer to the challenge of super-fast Chernoff approximation formulated by Remizov in Problem 30 of arXiv:2410.00416.

Nonlinear coupled dispersive models arise in numerous physical settings, including quantum electrodynamics, nonlinear optics, and plasma dynamics. In the nonrelativistic limit, such systems exhibit rapid temporal oscillations of wavelength $\mathcal{O}(\varepsilon^{2})$, controlled by a dimensionless parameter $0<\varepsilon\ll 1$. Moreover, owing to the Yukawa-coupled interaction terms, numerical schemes that simultaneously achieve high-order uniform accuracy and approximately preserve the mass and the energy are essential for reliable simulations. In this work, we propose a class of fully discrete schemes for nonlinear coupled dispersive models in the nonrelativistic regime that achieve higher-order uniform accuracy and approximately preserve both invariants over long times. Spatial discretization is based on the Fourier pseudospectral method, and time advancement employs two-scale exponential Runge-Kutta methods. Numerical experiments on the highly oscillatory Klein-Gordon-Dirac and Klein-Gordon-Schrödinger equations in one, two, and three dimensions confirm that the resulting schemes attain the theoretically predicted uniform accuracy and the long-time near-conservation of both invariants.

We propose a conservative local neural operator (CLNO) for accurate and generalizable autoregressive prediction of compressible flows. Building upon the locality and translation invariance of the local neural operator (LNO), CLNO introduces a conservative evolution framework based on local subdomains of fixed physical size, separating local mean-state evolution from fine-scale reconstruction. The flow is encoded into subdomain-integrated conservative states and fine-scale features, while a shared local network predicts interface fluxes for finite-volume updates that enforce discrete conservation. A convolutional reconstruction network recovers fine-scale structures using zero-mean basis functions, preserving the updated local means and providing information for subsequent evolution. To accommodate different physical parameters, Reynolds–Mach number conditioning modulates the encoding, flux-prediction, and reconstruction networks, enabling a single model to learn parameter-dependent dynamics and generalize to unseen parameter combinations without retraining. CLNO supports prediction across spatial resolutions, computational domain sizes, and boundary conditions. Numerical experiments demonstrate improved generalization and competitive or higher accuracy compared with LNO, including order-of-magnitude error reductions in representative cases, at approximately three to five times the inference cost on periodic benchmarks. These results demonstrate accurate, generalizable, and physically consistent flow prediction through conservative local evolution and fine-scale reconstruction.

We develop explicit second-order particle-in-cell (PIC) methods for the Vlasov–Ampère and Vlasov–Maxwell systems using the scalar auxiliary variable (SAV) approach. Two schemes are proposed for each system: a predictor-based scheme and a symmetric splitting scheme using a scaled electric field. The methods exactly preserve the discrete Gauss law and a modified total energy. Numerical experiments for electrostatic and electromagnetic problems demonstrate the expected dynamics as well as the conservation properties of both invariants.

Horizontally explicit, vertically implicit (HEVI) time integration removes the acoustic time-step restriction imposed by fine vertical grids in nonhydrostatic atmospheric models, but requires fast solvers for the resulting vertical linear systems, especially on GPUs. We present GPU-accelerated nested Schur-complement solvers for these systems for a split-form discontinuous Galerkin spectral element method on Legendre-Gauss-Lobatto (LGL) nodes applied to the compressible Euler equations with gravity. In time, we use a Rosenbrock-W method, a linearly implicit Runge-Kutta scheme that retains its order and stability with approximate Jacobians. This allows a sparse vertical Jacobian that captures vertical acoustic waves and gravity on curvilinear meshes with orography and decouples into independent column systems. The solvers are nested: a global Schur complement (static condensation) exploits the LGL sparsity to condense each column onto a few element-edge unknowns, and a local Schur complement eliminates all variables but the vertical momentum within each element, exploiting that the pressure depends only on the potential-temperature density. We propose and analyze two condensation strategies, which lead to banded and block-tridiagonal reduced systems. The resulting HEVI schemes conserve total mass and potential-temperature density to round-off. The solvers are implemented for CPUs and GPUs in TrixiAtmo.jl, the atmospheric extension of Trixi.jl. We verify high-order convergence with linear inertia-gravity waves in a 2D channel and demonstrate robustness on the sphere with the baroclinic instability, a mountain-induced baroclinic instability, and the DCMIP-2025 gap-flow and vortex-shedding cases. For the baroclinic instability without orography, the scheme achieves 18.8 and 2.9 simulated years per day at 100 and 50 km resolution on one NVIDIA GH200 GPU in double precision.

In this paper, we develop a greedy double-block nonlinear Bregman-Kaczmarz method to efficiently find sparse solutions of nonlinear equations. The theoretical convergence rate of the proposed method is established under certain assumptions. Numerical experiments illustrate the competitive performance and superiority of the proposed method.

High-order direct arbitrary Lagrangian–Eulerian discontinuous Galerkin (ALE-DG) methods can suffer severe order reduction under nonsmooth mesh motion despite satisfying the geometric conservation law (GCL). To address this limitation, we propose a new ALE-DG framework applicable to moving triangular and quadrilateral meshes. This framework successfully overcomes the uniform Lipschitz constraint inherent in existing direct ALE-DG methods, retaining high-order accuracy even when the mesh-velocity gradients are unbounded. This capability is driven by transport polynomial exactness (TPE($k$)), which is a structural property stronger than the GCL. TPE($k$) ensures the exact transport of polynomials of degree at most $k$ under arbitrary mesh motion. To incorporate TPE($k$) into a direct ALE-DG method, we introduce a coefficient-evolution (CE) formulation. Reconciling strict local conservation, an indispensable requirement for accurate shock capturing, with the preservation of the TPE(\(k\)) property constitutes the principal technical challenge. We resolve this by evolving the coefficients using the third-order strong-stability-preserving Runge–Kutta method, paired with a conservative correction executed via Simpson's rule integration. Crucially, this correction strictly vanishes for exactly transported quadratic polynomials, thereby ensuring that our third-order scheme rigorously satisfies the TPE($2$) property. Numerical experiments demonstrate that the proposed scheme retains third-order accuracy when the discrete mesh-velocity Lipschitz estimate grows as $\mathcal O(h^{-1})$, whereas a representative existing direct ALE-DG scheme exhibits order reduction under the same mesh motion. Tests with discontinuous solutions further demonstrate robust and high-resolution performance under random and Lagrangian mesh motions.

We consider the recovery of functions of the weighted Wiener classes $W_s^\mu$ from perturbed Fourier–Hermite coefficients in the metric of $L_\infty(\mathbb R)$. For the summation methods of the class $\mathcal S^\theta$ we find error estimates on $W_s^\mu$. The estimates take into account the behaviour of the Hermite functions both in the oscillatory region and near the turning point: the maxima of functions with different indices are attained at different points, and the supremum of a sum over a block of indices is substantially smaller than the sum of the suprema. Owing to this, in particular, we establish a common lower bound for all finite diagonal methods; for $s\ne4/3$ and $p\ne4/3$ it matches the upper bound and gives the exact order of the error in this class, while at the critical points the gap between the bounds does not exceed $(\log(e/\delta))^{1/4}$ up to a constant. If, in addition, $s>4/3$ or $p<4/3$, then this order is also optimal among all recovery operators that use the whole sequence of data.

Projected Walk on Spheres (PWoS) solves surface PDEs through ambient sampling and closest-point projection, but tightly folded embeddings can force short steps or ambiguous projections despite simple intrinsic geometry. We introduce Walk-on-Ambient-Patch (WoAP), which treats the embedding as a solver variable. Edge-isometry constraints and tangent-point repulsion separate intrinsically distant regions, optionally in higher codimension; patch states track residual projection branches without introducing PDE boundaries. We analyze the resulting geometry–sampling tradeoff: additional dimensions provide separation freedom but reduce uniform tangent progress, while strong tangent-biased proposals incur large importance-weight moments. Folded-surface experiments demonstrate reduced projection conflicts, a separate hard-constraint pilot verifies metric preservation, and an internal fixed-RMSE comparison favors higher dimensions at stricter accuracy targets. Together, these results motivate choosing the embedding and projection state jointly with the stochastic solver.

Nonlinear Constrained Optimisation (NCO) requires simultaneous satisfaction of the Karush-Kuhn-Tucker (KKT) conditions. Existing methods enforce these conditions through projections, penalties, barriers, or sequential approximations, potentially introducing computational complexity and convergence difficulties. For nonconvex problems, a KKT point need not be a local optimum. This paper reformulates NCO as stabilisation of a unified KKT sliding residual manifold and develops SM-KKT, a projection-free solver. Stationarity, equality feasibility, and regularised Fischer-Burmeister complementarity are integrated into a unified residual system, with complementarity encoding inequality feasibility, dual feasibility, and complementary slackness. The manifold's zero set is shown to be equivalent to the KKT conditions. The required Jacobian saddle-point inertia is characterised and shown, under stated assumptions, to be equivalent to the simultaneous satisfaction of the Linear Independence Constraint Qualification (LICQ) and Second-Order Sufficient Condition (SOSC). This characterisation establishes local stability of the structured correction algorithm under these conditions and distinguishes nondegenerate locally optimal KKT solutions from points failing second-order optimality. SM-KKT computes structured primal–dual corrections and applies a discrete sliding reaching law to drive the KKT residual into a prescribed boundary layer in finite time. Four nonlinear benchmarks show recovery of reference constrained optima with KKT residuals below $10^{-6}$, earlier boundary-layer reaching, and favourable KKT accuracy compared with projected primal-dual, augmented-Lagrangian, quadratic-penalty, and log-barrier baselines.

We derive polynomial bubble duality in barycentric coordinates and construct explicit nodal bases for the first family of finite element differential forms. Logarithmic coordinates give an elementary proof of the Whitney–bubble decomposition. The Hessian of negative entropy defines a metric whose normalized Hodge star maps bubble generators to complementary-degree polynomial forms. This gives a direct geometric description of the known dualities and proves nondegeneracy of their integral pairings. Combining Whitney forms with shifted Lagrange polynomials then gives bases dual to tangential point evaluations in every dimension and order. The construction extends to scalar nodes that are unisolvent on every face. The interpolant is independent of the roots and commutes with face restriction.

Stochastic gradient descent (SGD) and randomized block coordinate descent (RBCD) belong to two fundamentally different classes of stochastic iterative regularization methods, and are promising and popular for solving large-scale inverse problems, due to their excellent scalability with respect to data size and signal size, respectively. In this work, we propose a doubly stochastic block gradient descent (DSBGD) method with early stopping, by integrating the updating strategies of both SGD and RBCD, for solving linear inverse problems in infinite-dimensional Hilbert spaces. At each iteration, the method randomly selects one data pair to obtain a stochastic gradient estimate, and then uses it to update one randomly selected block coordinate. We prove the regularizing property of DSBGD under a priori stopping rules and then derive its convergence rates with the additional source condition and a condition on the forward operator. To the best of our knowledge, these are the first convergence rate results of randomized coordinate descent-type iterative methods for inverse problems from the perspective of regularization theory. We also prove that the DSBGD iteration can be stopped by the classical discrepancy principle in finitely many steps with high probability, and that the mean squared error decreases monotonically before stopping. Several numerical experiments are provided to complement the theoretical analysis.

In this paper, we propose and analyse a novel explicit fourth-order symmetric low-regularity trigonometric integrator for the semilinear Klein-Gordon equation with non-smooth solution in the $d$-dimensional space, where $d=1,2,3$. The integrator is constructed based on the full use of Duhamel's formula, the employment of a twisted function tailored for trigonometric integrals and the symmetrization procedure of trigonometric integrators. Robust error analysis is conducted, demonstrating that the proposed scheme achieves fourth-order accuracy in the energy space under a weak regularity requirement in $H^{2+\max(\mu,1)}(\mathbb{T}^d)\times H^{1+\max(\mu,1)}(\mathbb{T}^d)$ with $\mu> \frac{d}{2}$. A numerical experiment shows that the proposed low-regularity integrator is much more accurate than some well-known exponential integrators for approximating the Klein-Gordon equation with non-smooth solutions.

It is well known that the matrix structure affects the convergence of iterative methods, particularly those utilizing Krylov subspaces. The matrix structure, determined by its graph, dictates the propagation order of updates to the vector components that are recursively computed during the conjugate gradient method iterations. Along with the propagation of these component updates, round-off errors also propagate. We present a model for the propagation and accumulation of round-off errors and demonstrate that the larger the values of topological characteristics of the matrix graph, such as its diameter, effective diameter, and pseudodiameter, the slower the round-off error accumulates during the method's iterations, and the longer the propagation of round-off errors remains local when computing the components of the search direction, residual, and approximate solution vectors. Round-off errors arising during the computation of these components affect only those components whose corresponding vertices in the matrix graph are at a distance not exceeding the current iteration number. As the iteration number increases, the propagation of round-off errors between the vertices of the matrix graph becomes global, leading to the loss of orthogonality of residuals and conjugacy of search directions. In the context of mixed-precision computations, this implies that for a matrix with a larger graph diameter (or, more practically, effective diameter), a greater number of method iterations can be performed in a lower-precision mode before the onset of convergence stagnation.

We develop an adaptive Galerkin neural-network method for computing minimum-$L^2$-norm interior controls for wave equations with spatially varying coefficients. Within the Hilbert Uniqueness Method, adjoint waves are approximated directly in space-time by least-squares Galerkin projections onto trial spaces enriched by neural-network basis functions selected through residual-based training. Linear Galerkin solves determine the adjoint approximations, and a separate linear system reconstructs the control. The least-squares objective provides an indicator for adaptive enrichment and stopping; the experiments use a Ritz-based terminal-data version whose continuous-error interpretation additionally requires an eigenpair consistency estimate. The main difficulty is that adjoint errors perturb this reconstruction system. Assuming observability of the continuous equation, we derive sufficient stability conditions involving a worst-case factor $\sqrt{n}$ multiplying the maximal adjoint error, where $n$ is the number of retained terminal modes. A fully discrete error estimate separates terminal-space truncation, adjoint approximation, and assembly errors, yielding convergence to the minimum-norm control when the truncation error, the $\sqrt{n}$-weighted adjoint error, and the assembly errors tend to zero as $n$ grows. For spectral terminal bases, additional regularity and uniform contraction assumptions yield quantitative convergence rates and show that $O(\log n)$ accepted enrichment steps per mode suffice to balance the neural and truncation errors. Numerical experiments on a nonconvex Cassini oval, a square containing a low-speed lens, and a three-dimensional ball with a graded coefficient illustrate the method.

We study a local graph-based interpolation framework for data defined on graph vertices embedded in a metric space. The resulting interpolants are constructed from variational graph splines defined on small neighborhoods, yielding continuous-domain extensions of the graph data. We analyze their behavior through energy estimates, Lipschitz continuity, and stability properties, and compare them with $k$-Nearest Neighbor ($k$NN) interpolation as a natural baseline due to its simplicity and widespread use. Numerical experiments on the manifold of symmetric positive definite matrices illustrate the relative performance of the proposed local interpolation schemes and $k$NN interpolation across varying levels of smoothness and noise.

In deep hierarchies of adaptive mesh refinement, the update work that subcycling avoids is set mainly by the dimension of what is refined, not of the space. Subcycling (local time stepping) does less work than global stepping, advancing each level with the finest step, provided its extra coordination cost stays below an overhead budget $\beta=1-1/\Sid$ as a fraction of the global-stepping work. Here $\Sid$ is the ideal speedup that the cell counts give for equal update costs. When each level follows a feature of box-counting dimension $\delta$ and the counts grow geometrically, the budget of a code that advances every cell tends to $(r-1)/(r^{\delta+1}-1)$ for refinement ratio r. A point leaves a budget near one, while a surface in three dimensions leaves 1/7 at r=2 and 1/21 at r=4, values implied by the planar estimates of Dursi and Zingale. Published comparisons in discontinuous Galerkin and octree codes with local time stepping imply overheads of 0.1 to 5 percent of the work, below 1/7 and comparable with 1/21. In paired runs of one AMReX code on one core, a cost model fixed in advance and fitted on points and filled regions predicted the speedups of curves and surfaces with a median error of 7 percent (22 percent on the saving); post hoc, the update counts alone did as well except for points. The overhead, as defined here, was negative in all 17 configurations, but when computed from speedups like the published values, it ranged from -0.11 to 0.19, and no configuration came near break-even. When an update costs the same under both schedules, one profiled run can show that subcycling pays.