Papers

The low-rank alternating direction implicit (ADI) method is an efficient numerical technique for solving several types of large-scale matrix equations that admit low-rank solutions. The discrete-time algebraic Riccati equation (DARE) is an important matrix equation with applications in state estimation, controller design, and filter design. In the literature, the low-rank Cholesky factor ADI method for Stein equations has been used within Newton iterations to solve large-scale DAREs. However, no dedicated low-rank ADI solver is available for such DAREs. To address this gap, this paper presents a low-rank ADI solver for large-scale DAREs. We also propose an efficient approach to generate ADI shifts automatically, which makes the proposed solver fully autonomous for solving DAREs. The effectiveness of the proposed solver is compared with MATLAB's idare on a moderate-order problem. Efficiency and accuracy are further demonstrated on large-scale DAREs of order $10^6$. Numerical results confirm that the solver is efficient, accurate, and fully autonomous.

This paper presents a convergence analysis for a newly developed nonlinear model reduction method: parametric probabilistic manifold decomposition (PPMD)~\cite{guo2026parametric}. In addition, existing analyzes of nonlinear reduced order models typically treat subspace reduction, manifold representation, regression, and nonlinear reconstruction as separate components and often remain at the level of discrete state vectors. To the best of our knowledge, no theory tracks the complete error propagation in a data-dependent model whose basis, residual geometry, spectral coordinates, parameter maps, and lifting operator are all learned from the same numerical solution data. We develop a coupled perturbation analysis for the entire PPMD procedure. A trajectory geometry induced by the spatial discretization and temporal quadrature connects discrete trajectory vectors isometrically with the corresponding PDE norm. Population spectral objects are introduced to align the empirical residual coordinates and derive a uniform coordinate error estimate, whose propagation through the Hilbert-valued kernel lifting estimator is then quantified. Combining these results with the full order discretization error, weighted low-rank approximation, parameter regression, and residual representation defect yields deterministic and high-probability trajectory error bounds and consistency in probability in the continuous PDE trajectory space. The theory identifies how the principal errors interact and which components limit the accuracy of the nonlinear reduced order model.

This paper presents a novel framework that integrates learnable graph convolutions with the geometric multigrid method for solving partial differential equations (PDEs). The discretization of PDEs is first represented as a graph structure, enabling the application of graph convolutions to enhance the multigrid performance. By incorporating graph convolutions into the multigrid components such as smoothing and inter-grid transfer operators, we develop a learnable multigrid that can adaptively optimize its performance based on the underlying problem characteristics. In this framework, the graph convolutions are embedded directly within the multigrid cycle, effectively transforming the entire multigrid solver into a specialized neural network architecture, rather than combining a classical solver with surrogate models. The learnable multigrid framework is lightweight in terms of parameter count and requires minimal training effort to achieve good performance. Numerical experiments demonstrate the effectiveness of the proposed approach in solving some challenging problems, showing improved convergence rates compared to traditional multigrid methods. The generalizability of the learned parameters across different problem settings, including varying source terms, coefficients, geometries, and mesh sizes, is also investigated with proper weight-sharing and transfer-learning strategies.

This paper proposes gappy probabilistic manifold decomposition (Gappy PMD), a nonlinear method for reconstructing high-dimensional fields from extremely sparse measurements. Gappy PMD reconstructs the field on the nonlinear manifold learned by probabilistic manifold decomposition (PMD). We further propose a differentiable point selection method for reduced-order model (ROM)-based field reconstruction (DPS). Using differentiable meshless interpolation within the ROM-based reconstruction framework, DPS makes the full-field reconstruction error differentiable with respect to the sampling locations and directly optimizes these locations. In addition, a theoretical error analysis for Gappy PMD is also given. It splits the squared reconstruction error into two orthogonal parts: one normal to the reconstruction manifold and the other induced by sparse sampling and observation noise. Under a stability condition on the sampling operator, this error vanishes with the PMD approximation error and the noise. The Gappy PMD is evaluated on three numerical test cases: flow past a cylinder, lid-driven cavity flow, and backward-facing step flow. For the same reduced dimension and sampling points, Gappy PMD attains mean relative $L^2$ errors one to two orders of magnitude below Gappy POD. Optimizing the sampling points with DPS further improves reconstruction accuracy and robustness.

Let $$ A=\begin{bmatrix} H_1 & E^* \\ E & H_2 \end{bmatrix} \quad\text{and}\quad \widetilde A=\begin{bmatrix} H_1 & 0 \\ 0 & H_2 \end{bmatrix} $$ be partitioned Hermitian matrices, with $\widetilde A$ obtained by dropping the off-diagonal blocks of $A$. Let $\eta$ be the minimum distance between the spectra of $H_1$ and $H_2$; these spectra may interlace. Define, for $\delta,\epsilon\ge0$, $$ \phi(\delta,\epsilon)= \begin{cases} 2\epsilon/(\delta+\sqrt{\delta^2+4\epsilon^2}), &\quad\mbox{if $(\delta,\epsilon)\ne (0,0)$}, 1, &\quad\mbox{if $(\delta,\epsilon) = (0,0)$}, \end{cases} $$ and let $\widetilde E=A-\widetilde A$ and $\epsilon_2=\|E\|_2=\|\widetilde E\|_2$, the matrix spectral norm. Li and Li [Linear Algebra Appl., 395:183–190, 2005] established a sharp spectral-norm bound on the changes in the eigenvalues of $A$: $$ \big\|\mathrm{diag}\big(\pmb{\lambda}(A)-\pmb{\lambda}(\widetilde A)\big)\big\|_2 \le \phi(\eta,\epsilon_2)\,\|E\|_2, $$ where $\pmb{\lambda}(A)$ is the vector whose components are the eigenvalues of $A$ in descending order and similarly for $\pmb{\lambda}(\widetilde A)$. In this paper, this bound is extended to every unitarily invariant norm $\|\cdot\|_\mathrm{ui}$: $$ \big\|\mathrm{diag}\big(\pmb{\lambda}(A)-\pmb{\lambda}(\widetilde A)\big)\big\|_\mathrm{ui} \le \phi(\eta,\epsilon_2)\,\|A-\widetilde A\|_\mathrm{ui}. $$ A $2\times2$ example attains equality for every such norm. An analogous bound for singular values is also established when the off-diagonal blocks of a rectangular matrix are dropped.

Many open systems are observed only through aggregate snapshots, while the rule that moves objects between states stays hidden. This paper studies the inverse problem of recovering the transfer operator of a multicluster open system from two adjacent snapshots. The consistent operators fill an affine manifold: the data fix the aggregate flows and leave the internal routing free, so every recovery method commits to a prior that selects one operator. Tikhonov regularisation with a zero reference gives the uniform outflow on each source simplex, which underestimates retention in an inertial system. We propose least action on the operator trajectory, under which the rule changes as little as the data permit. In discrete time this is exact as a probability statement: the action prices a trajectory by its likelihood under a Gaussian reference drawn toward full retention, so least action selects the most probable trajectory consistent with the balance. A causal construction keeps the balance exact by moving the operator only inside the null space of the constraints, and we relate it to optimal transport, to Perron-Frobenius coarse-graining, and to the Schrödinger bridge. Two instantiations carry it: the Russian labour market and a radiation-damped relativistic electron ensemble. The construction holds the balance to machine precision where temporal smoothing corrupts it by twelve percent on the labour series and sixty percent in the physical one. The spectral gap is non-identifiable from aggregate snapshots but available from micro-trajectories, where Ulam's method reproduces the Landau-Lifshitz rate to within six percent. On the short labour series no method beats the naive forecast in point accuracy at conventional significance levels. The value is structural: an interpretable operator, an exactly preserved balance, and an explicit information limit of aggregate observation.

We prove, to our knowledge, the first unconditional cumulative closure estimates for the Arnold–Mukherjee–Pouly (AMP) refinement algorithm and the original face-marked tetrahedral algorithm of Bänsch on arbitrary conforming initial tetrahedral meshes. Let $T_0,\ldots,T_L$ be an adaptive mesh sequence generated by either algorithm, with $M_\ell\subseteq T_\ell$ denoting the marking set at step $\ell$. Then $$\#T_L-\#T_0\leq C_{\mathrm{clos}}(T_0)\sum_{\ell=0}^{L-1}\#M_\ell.$$ The proof is intrinsic to the physical three-dimensional mesh and requires neither an initial compatibility condition nor a higher-dimensional embedding. It organizes conformity refinements into a causal forest and combines a uniform horizontal-propagation estimate with a weighted packing argument to obtain an explicit closure constant. For the original Bänsch algorithm, every history-dependent resolution of the initial two-edge ambiguity is represented by one of finitely many AMP histories. The estimate therefore holds uniformly for arbitrary deterministic or nondeterministic choices. This resolves a long-standing complexity question for the Bänsch–AMP family.

We present a numerical formulation for the resolution of finite systems of interacting edge dislocations and wedge disclinations. The approach is based on a finite element discretization of a fourth-order elliptic boundary value problem arising from the mechanical equilibrium equations of plane-strain linear elasticity in the presence of kinematic incompatibilities. The numerical implementation follows a continuous interior penalty discontinuous Galerkin framework and relies on the solution of a finite set of local cell problems. Our method is validated against analytical benchmarks and used to explore several non-trivial dislocation-disclination configurations, illustrating how translational and rotational incompatibilities shape the stress state of the body.

We establish a uniform inf-sup stability estimate for the Scott-Vogelius finite element spaces on uniform Freudenthal tetrahedralizations of the unit cube for polynomial degrees k >= 4. This result completely settles the first conjecture of Farrell, Mitchell, and Scott for the critical degrees k = 4 and k = 5, complementing the known stability range for higher polynomial degrees. The main mathematical difficulties stem from the complex topological compatibility required at the singular vertices and the corresponding mean-value constraints across adjacent elements. We tackle these challenges by developing a unified barycentric skeleton-bubble calculus that explicitly constructs vertex jets, edge modes, and face transfers to globally route element means. The accompanying exact computations independently verify these finite-dimensional identities and provide reproducibility data.

Solving ill-conditioned linear systems needs LU decomposition in precisions beyond binary64. The standard approaches – GMP/MPFR, or multi-component arithmetic such as double-double – leave every scalar multiply-add inside the O(n^3) update in multiple precision, and so cannot exploit the low-precision matrix engines that dominate current hardware. We build a blocked, partially pivoted LU decomposition on Ozaki scheme II, which replaces the multiple-precision GEMM by exact integer modular products followed by an explicit CRT reconstruction, and implement it for multi-component (DD/TD/QD) and arbitrary precision, on CPUs and GPUs. Two contributions make this practical: a direct conversion between the non-overlapping expansion format and the internal fixed-point representation, which removes the MPFR round trip and is worth a factor of 2.4–4.7; and new FP16, FP8 and binary64 GPU back-ends, the FP8 one using a balanced base-17 two-digit encoding that keeps the full 362.8-bit CRT capacity of INT8. On an Arm/GB10 and an x86/H100 the fastest back-end changes with the machine: INT8 wins on GB10, whereas on H100 binary64 is fastest for QD, beating both the native implementation and INT8. The essential point of Ozaki scheme II is thus not to use a low-precision engine but to choose the format maximising bits-per-modulus times engine throughput. For the Lotkin matrix at p ~ 1.2 log2 cond(A) we reach relative errors of 10^-645 at n=2048, up to 2.84x faster than a fully OpenMP-parallel multiple-precision LU. We also show by measurement that the O(p) advantage of scheme II over scheme I applies only to the GEMM term: modular reduction and CRT grow as O(p^2) and dominate the runtime at the matrix sizes considered here. The implementation is released as open source.

This paper presents several parallel versions of Jacobi's method for the symmetric eigenvalue problem and the SVD. A continuation of [Demmel, Luo, Schneider, & Wang 2025], we develop parallel Jacobi algorithms whose arithmetic cost is optimal and whose bandwidth or latency can match the corresponding lower bounds of parallel matrix multiplication. Our focus is a standard distributed-memory setting with variable processor layouts, including both 2D and 2.5D processor grids. In the 2D case, we demonstrate that a standard implementation of parallel Jacobi achieves a perfect speedup in arithmetic cost – i.e., complexity $O(n^3/P)$ when done with $P$ processors – while hitting the 2D matrix-multiplication lower bound for bandwidth and (nearly) the lower bound for latency. By employing a 2.5D processor grid and leveraging 2.5D matrix multiplication, equivalently by increasing the memory per processor, we demonstrate that parallel Jacobi can achieve even lower bandwidth/latency, though we also prove that these costs cannot simultaneously match the best-known bounds for parallel matrix multiplication in any Jacobi algorithm. Finally, we extend our results to one-sided Jacobi SVD.

Low-rank methods can reduce the memory and computational requirements of deep neural network training in approaches such as GaLore. Randomized range finders offer an attractive alternative to singular value decompositions, particularly when the required rank is determined adaptively from a prescribed approximation tolerance. We introduce a blocked adaptive randomized range finder based on an implicit Householder QR decomposition and an optimized hybrid CPU–GPU implementation. The proposed method avoids explicit reorthogonalization. Numerical experiments show that it preserves orthogonality and approximation accuracy in regimes where block Gram–Schmidt without reorthogonalization becomes unstable. The blocked formulation exposes matrix–matrix operations and enables overlap of CPU panel factorization with GPU updates. On an NVIDIA GH200, the overlapped implementation reduces the runtime for the largest tested matrix from 9.91 seconds on the CPU to 0.407 seconds. The method provides a stable and efficient building block for low-rank approximation on heterogeneous systems with applications in computational science and engineering.

Cubature rules, which approximate definite integrals as a linear combination of a set of function values, are ubiquitous and necessary for computational methods in the physical sciences. A superconvergence result for cubature rules on polytopal domains in two and three dimensions is developed, whereby a rule that is exact for all bivariate polynomials of a fixed even degree realizes an extra order of convergence under a decrease in the spacing between nodes.

In Nyström-collocation discretizations of the electric field integral equation (EFIE), the surface divergence acts on surface densities that may be discontinuous across patch boundaries, which degrades accuracy and convergence. We show that this not only affects the EFIE but every formulation in which the operator occurs, either in the equation itself or in the scattered field computation, and propose a high-order-accurate continuity-enforcing scheme for smooth surfaces as a remedy for the direct and indirect EFIEs, magnetic field integral equations (MFIEs), and regularized combined field integral equations (CFIEs) alike. The scheme comprises two ingredients: i) We show how to discretize the equations via a Chebyshev-based Nyström scheme, which admits closed quadrature rules. ii) Since unknowns and test vectors are in terms of patch-local curvilinear bases, continuity is enforced by a change of basis: we construct sparse mapping matrices assembled solely from the curvilinear geometry description. In doing so, we restore the accuracy of the EFIE such that it can be combined with the MFIEs with equal weights to form CFIEs. Numerical studies for the scattering from canonical and realistic geometries show that all considered formulations individually and combined benefit from the continuity enforcement in terms of better conditioning, reduced iterations of an iterative solver, and several more digits of accuracy in the scattered fields, despite reducing the total number of unknowns.

Separable nonnegative matrix factorization (SNMF) has been widely used for low-rank representation and clustering of nonnegative data, owing to its ability to produce part-based and interpretable decompositions. In particular, SNMF is closely related to graph clustering and community detection. To enhance sparsity and identifiability of the learned factors, we propose an $\ell_1^p/\ell_2$-regularized SNMF model based on a powered ratio-of-norms regularizer. The resulting formulation is nonconvex and nonsmooth, which poses significant challenges for optimization. To address this, we develop efficient algorithms based on the difference-of-convex function algorithm (DCA) and the alternating direction method of multipliers (ADMM). The proposed methods decompose the original problem into tractable subproblems, leveraging closed-form proximal operators associated with the powered norm terms. We establish descent and limiting criticality properties for the DCA scheme and convergence under standard assumptions for the ADMM scheme. Extensive numerical experiments on synthetic datasets and hand gesture classification tasks demonstrate that the proposed approach achieves competitive or improved performance in anchor identification and classification accuracy compared with existing SNMF methods, while maintaining competitive computational efficiency.

We introduce a fully discrete numerical scheme for semilinear SPDEs with additive space-time white noise that overcomes the previous order barrier for the spatial convergence rate. The scheme, which we refer to as the doubly accelerated exponential Euler scheme, achieves a strong convergence rate of $M^{-1+\epsilon}$ in time and $N^{-3/2+\epsilon}$ in space for any $\epsilon>0$, where $M^{-1}$ and $N^{-1}$ are the temporal, respectively the spatial, meshsizes. This substantially improves the standard spatial error bounds of order $N^{-1/2}$ in the literature. Numerical simulations support the findings.

We propose QGPINNs, a physics-informed neural network framework developed in PyTorch for the numerical solution of nonlocal differential equations on quantum graphs. The framework is designed as a general computational implementation in which the solution on each edge of the graph is approximated by a neural network, while a unified graph-based loss function enforces the governing equations together with initial, boundary, and vertex transmission conditions. In particular, the formulation incorporates standard continuity and Kirchhoff-Neumann vertex conditions and Dirichlet boundary conditions into the learning process to couple the local edge-wise neural approximations into a global solution on the graph. The framework is developed for two representative classes of nonlinear models: multi-order fractional elliptic problems and time-fractional evolution equations on quantum graphs. To improve accuracy and training stability, QGPINNs integrates several graph-adapted learning strategies, including soft and hard constraint enforcement, dynamic loss balancing, Fourier feature embeddings, and a learnable singularity-capturing feature for weakly singular solutions arising in the considered problems. The framework also extends naturally to inverse problems, including the identification of the orders of fractional operators and physical parameters from noisy observational data. We validate the accuracy, computational efficiency, and physical consistency of the proposed framework through numerical experiments on benchmark graph structures and real-world networks, including the IEEE 14-bus system and an open-channel agricultural drainage network.

We propose an approach based on Riemannian optimization to compute a nearest normal matrix to a given one. The problem can be formulated as the minimization of a smooth function either on the manifold $U(n)$ of unitary matrices of size n or on the flag manifold $U (n)/U (1)^n$. The flag manifold is particularly suitable for theoretical analysis; we characterize the global maximum of the objective function and prove that, for generic inputs, its local minimizers are finitely many and isolated; in turn, this implies the original nearest normal matrix problem generically has finitely many local minimizers, all with distinct eigenvalues. We also develop a Riemannian trust-region method that improves substantially on classical algorithms and can handle considerably larger matrices, as well as a variant for computing the nearest real normal matrix. The paper is complemented by extensive numerical experiments.

Neural operators provide fast surrogate models for approximating operators between function spaces, but their predictions often lack uncertainty quantification. We develop a split conformal framework to guarantee that a calibrated pointwise band around the neural operator output contains the true solution on at least a $1-\gamma$ fraction of the evaluation domain, with probability at least $1-\alpha$ over test and calibration inputs, where $\alpha,\gamma\in(0,1)$. Our method reduces a normalized residual field to its spatial $(1-\gamma)$-quantile and computes a scaling factor using a held-out calibration dataset. We prove marginal coverage guarantees for measurable residual fields defined on arbitrary probability spaces, covering both continuum domains and fixed discretizations. Under mild assumptions on the data distribution, we show that the coverage conditional on the calibration set follows a Beta distribution, which we verify with numerical experiments on Darcy flow and Navier–Stokes equations, where our calibration yields bands consistently tighter than existing corrections while retaining the target coverage.

This work investigates a time-dependent source identification problem for linear parabolic equations subject to an integral constraint and Neumann boundary conditions in a domain of $\mathbb{R}^d$, $d\ge 1$. We establish well-posedness and higher regularity of the solution pair in parabolic Hölder spaces. A numerical algorithm based on a finite element discretization in space and an implicit time-stepping scheme is then developed for the reconstruction of the unknown source. The resulting discrete inverse problem involves a well-conditioned operator. We show that identity Tikhonov regularization provides only uniform, nonselective shrinkage in this setting, whereas Tikhonov regularization with derivative-based penalties, analyzed through the generalized singular value decomposition, provides an effective denoising strategy. The regularization parameter is selected using the Morozov discrepancy principle. Numerical errors are evaluated using full parabolic Hölder norms, which provide a more comprehensive assessment of the reconstruction by incorporating errors in the solution, its derivatives, and the associated Hölder seminorms. Numerical experiments for smooth and piecewise constant sources demonstrate accurate and robust reconstructions under increasing levels of noise.

Implicit trace estimation aims to approximate the trace of a matrix or linear operator accessible only through matrix-vector or operator-vector products. In the matrix setting, the Girard-Hutchinson estimator typically requires $\mathcal{O}(\varepsilon^{-2})$ products to achieve accuracy $\varepsilon$, while the variance-reduced Hutch++ algorithm reduces this sample complexity to $\mathcal{O}(\varepsilon^{-1})$ for positive semidefinite matrices. We develop infinite-dimensional analogues of these estimators for positive trace-class operators on separable Hilbert spaces. The idealized estimators use Gaussian random elements whose covariance is determined by the target operator, leading to unbiased operator versions of Girard-Hutchinson and Hutch++. We prove high-probability error bounds analogous to the finite-dimensional matrix results; in particular, idealized infHutch++ achieves $\mathcal{O}(\varepsilon^{-1})$ sample complexity. For practical computation, we introduce truncated implementations that restrict the random samples to finite-dimensional subspaces; for fixed sample budget, we show that truncated infHutch++ converges in distribution to its idealized counterpart as the truncation dimension tends to infinity. Numerical experiments with integral operators, density-of-states approximations, and spectral filtering for a radial Dirac operator show that these truncated estimators can achieve accuracy comparable to the ContHutch++ algorithm by Zvonek, Horning & Townsend while using lower-degree function representations and smaller internal discretizations in chebfun.

We present a novel unfitted computational framework for simulating fully nonlinear potential flow-based water waves. Focusing on wave propagation, we describe the core methodology, which involves a high-order polynomial-corrected shifted-boundary approximation on unfitted spectral elements. This approach allows for the simulation of a curved, highly time-dependent (moving and deforming) free surface affected by bathymetric changes. All that on a simple Cartesian background mesh without re-meshing or further approximations of non-affine geometric features. In addition, we highlight the importance of proper gradient recovery using a polynomial-preserving technique to accurately capture the vertical free-surface velocity. Ultimately, the goal is to develop a high-order convergent numerical scheme capable of simulating highly nonlinear waves over long periods. This is achieved through an arbitrary-order finite-difference approximation of the free surface combined with added hyperviscosity for numerical stability. We present verification and validation test cases for wave propagation in both periodic and finite domains. Emphasis is placed on convergence studies, the justification for using high-order approximations, the importance of optimal gradient recovery, long-time simulation of highly nonlinear waves, and nonlinear wave interactions with both affine and curved bathymetry changes, as well as with vertical walls, for highly nonlinear stream function waves and high-amplitude solitons. The novel computational model offers a comprehensive, all-in-one framework for simulating ocean waves and their interactions with offshore structures. It provides significant geometric flexibility, enabling boundaries such as the free surface to move and deform over time without the need to re-mesh a boundary-fitted mesh.

This paper investigates an inverse source problem for a time-fractional diffusion-wave equation with a singular inverse-square potential. The source term is assumed to consist of a known temporal factor and an unknown spatial component, which is to be recovered from terminal-state measurements. The well-posedness and regularity of the forward problem are established within an appropriate energy framework by exploiting Hardy-type inequalities and the spectral properties of the associated singular elliptic operator. The terminal observation operator is then shown to be compact, and uniqueness of the spatial source is established under a suitable nondegeneracy condition on the temporal factor. To stabilize the resulting ill-posed inverse problem, a Tikhonov regularization approach is introduced. The gradient of the regularized functional is derived through an adjoint problem involving a right-sided fractional derivative, leading to an adjoint-based conjugate gradient method with an exact line search for the numerical reconstruction of the unknown source. Numerical experiments are conducted on both one-and two-dimensional spatial domains, using both exact and noisy terminal data, to demonstrate the effectiveness and stability of the proposed source reconstruction method.

In this paper, we propose a local discontinuous Galerkin (LDG) method for the modified Camassa-Holm (mCH) equation that contains cubic nonlinearity and high-order derivative terms. Energy conservation and stability are obtained using the conservative and dissipative fluxes, respectively. The error estimate of the semi-discrete scheme is also established. Numerical examples verify that our theoretical findings are sharp and the proposed LDG method is capable of capturing the peakon solution of the mCH equation.

This manuscript presents a performance evaluation of Fast Fourier Transform (FFT) implementations on emerging processors supporting the RISC-V Vector Extension (RVV 1.0). By introducing juFFTe, a light-weight high-performance library for discrete Fourier transforms, it is demonstrated how effective vectorization of performance-critical FFT kernels can be achieved on RVV-enabled hardware. Comprehensive benchmarks on three RVV 1.0-ready processors, the SiFive X280, the X100 core of the SpacemiT K3 and the C920v2 core of the Sophon SG2044, reveal substantial performance improvements of juFFTe (https://github.com/FZJ-JSC/juFFTe) over the widely used FFTW3 library. Although RVV-enabled platforms show promising results at this stage of development, a comparison with AMD's Zen 5 architecture indicates that RISC-V needs further maturing to reach the performance of established micro-architectures.

Retractions provide a computationally efficient alternative to the Riemannian exponential and logarithm maps for practical data-processing tasks on manifolds. In particular, second-order retractions with closed-form inverse are well-suited for interpolation problems on manifolds. On the Stiefel manifold of orthogonal frames, there are only two retractions of this type: the Cayley retraction, which is second-order accurate under the canonical metric, and the recently proposed polar-light retraction, which is second-order accurate under the Euclidean metric. In this paper, we study the properties of these maps in the context of interpolation on the Stiefel manifold. To obtain explicit interpolation error bounds, we examine the conditioning of the retraction maps and their inverses. We show that the retractions are well-conditioned, and we derive interpolation error bounds similar to those of classical Euclidean interpolation. The inverse retractions are not well-conditioned in general, and we discuss how data can be mapped via an isometric group action to ensure stable computations. As with all retractions on compact manifolds, the inverse canonical Cayley retraction and the inverse polar-light retraction exist only locally, and we construct normal neighborhoods around any point in which either the inverse Cayley retraction or the invese polar-light retraction are guaranteed to be computable. As an application of the retraction maps, we consider Hermite interpolation, where the objective is to reproduce both sampled function values and derivative information. A numerical example demonstrates that retraction-based interpolation is competitive with classical methods based on Riemannian normal coordinates.

Diffusion synthetic acceleration is most effective when its diffusion correction reflects the transport discretisation that generates the iteration error. We develop this principle for high-order upwind discontinuous Galerkin discretisations of discrete-ordinates transport on polytopic meshes. From the discrete transport sweep, we derive the exact scalar correction that removes the source-iteration scalar error in one step. We prove that the associated scalar response is positive and self-adjoint, obtain an exact expression for the source-iteration convergence factor, and quantify the additional damping produced by vacuum leakage. Using the exact correction as a reference, we construct a transport-matched modified interior penalty correction whose boundary terms are inherited directly from homogeneous vacuum inflow. In the optically thick regime, the resulting MIP form approximates the exact correction with relative error proportional to the effective cell Knudsen number. This gives a strict acceleration of source iteration, with bounds uniform in mesh size, polynomial degree, and element face count for admissible polytopic meshes. Numerical experiments on Cartesian and centroidal Voronoi meshes confirm the predicted convergence and correction-operator scaling.

Quantile randomized Kaczmarz (QRK) has proven to be an efficient solver for corrupted linear systems and has received much attention. It was recently shown by Cai et al. (SIAM J. Matrix Anal. Appl. 47(2):802-823, 2026) that using $O(\log T/\log(1/\beta))$ samples for computing the quantile is necessary and sufficient for QRK to converge linearly over $T$ iterations when solving linear systems with a $\beta$-fraction of arbitrary corruptions, as long as $\beta$ is small enough. However, it remains unclear how large the corruption level $\beta$ can be, and how to compute the required subsample size $D$ explicitly, without hidden constants. This paper studies streaming linear systems with Massart noise via QRK using an order-optimal batch size $D=O(\log T)$ in each update. The independence of samples from previous iterations in the streaming setting enables a sharper analysis, yielding explicit, computable bounds on both the tolerable corruption level and the required subsample size. In particular, we establish linear convergence for corruption levels of up to approximately 7%. We also discuss how the constants improve under oblivious noise.

In this work, we develop a two-level overlapping preconditioner for time-dependent problems in high-contrast multiscale media. We present a coarse-space construction based on multiscale methods, with emphasis on the relaxed nonlocal multicontinuum (NLMC) method. The NLMC space can be separated into components representing the high-permeability regions and the low-permeability background. We show that the complete NLMC space effectively preconditions the heterogeneous stiffness operator, whereas the high-permeability component alone captures the contrast-dependent global modes. For suitable small time-step sizes, the mass matrix controls the low-permeability contribution and only the high-permeability component in NLMC space is required for the global coarse correction in the two-level preconditioner. For general time-step sizes, performance can be maintained by using the full NLMC space or augmenting the high-permeability component with a standard multiscale space. The proposed coarse-space construction lowers the computational cost while preserving robustness with respect to coefficient contrast and fine-scale resolution. To further improve efficiency, the multiscale basis functions can be constructed by iteratively solving the relaxed energy-minimizing formulation. We demonstrate the robustness, efficiency and scalability of the proposed method through several numerical experiments.

We give a mathematical proof of the optimal convergence rates for the D2Q9 BGK lattice Boltzmann method with the half-way bounce-back rule for the incompressible Stokes equations in a flat channel. The convergence rates are second-order for the velocity and first-order for the pressure as the lattice spacing $h$ tends to zero, in agreement with formal analyses and numerical experiments, whereas the available rigorous convergence theorems only yield an $O(h^{1/2})$ bound for the velocity error. A key step in the proof is a decomposition of the leading boundary consistency error into macroscopic and kinetic components. These components are absorbed by suitably constructed Stokes and discrete Knudsen layer correctors, respectively. Incorporating these correctors into the prediction function used in previous rigorous analyses, we obtain a refined prediction function with consistency errors of sufficiently high order. Combined with the known weighted $L^2$-stability estimate, this gives the optimal convergence rates.