Papers

Mesh quality is very well studied and widely used to quantify a good mesh. In contrast, a systematic study of quality of meshfree point clouds is lacking. This gap makes it difficult to substantiate the common claim that generating a good-quality point cloud is easier than generating a good mesh. Various definitions of point cloud quality have been proposed, some of which have theoretical significance for proving convergence and error bounds, while others are used in computational studies. In this work, we compare and contrast existing point cloud quality metrics and introduce a few new ones. We conduct extensive numerical tests with a meshfree collocation method across a wide range of scenarios, including both elliptic and hyperbolic equations, 2D and 3D domains, and variations in parameters of the numerical method. Based on these tests, we assess which quality metrics best correlate with numerical error. Our findings reveal six metrics that consistently serve as reliable indicators of point cloud quality, while also demonstrating that several widely used metrics are poor predictors of accuracy.

The randomly pivoted Cholesky algorithm is one of the leading methods for computing a low-rank approximation to a large positive-semidefinite matrix. However, while it consistently achieves accuracy comparable to or better than competing methods of its type in experiments, its theoretical analysis lags somewhat behind other methods. This paper closes this gap, proving that randomly pivoted Cholesky produces an approximation with expected error within a $1+\varepsilon$ factor of the optimal rank-$r$ approximation in $\mathcal{O}(r/\varepsilon + r\sqrt{\log r})$ steps. This result nearly matches the optimal complexity $\Theta(r/\varepsilon)$ for any low-rank approximation method based on a partial Cholesky decomposition (also known as a column Nyström approximation). The paper also presents bounds on the randomly pivoted Cholesky trace and spectral-norm errors that hold with high probability. The mathematical argument is largely due to GPT 5.6-Sol (Pro), with some refinements by the author.

Turing completeness of a programming language or system characterizes its expressive power; and the strong Church-Turing hypo-/thesis refines such from qualitative to polynomial-time equivalence. Exact Real Computation (ERC) is a novel numerical programming language paradigm: for the imperative processing of continuous data as entities appearing as exact, i.e. devoid of rounding errors [doi:10.1007/978-3-662-44199-2_107]. ERC has been designed [doi:10.46298/lmcs-20(2:17)2024] as convenient and practical alternative, namely proven qualitatively equivalent, to the Turing machines originally underlying Computable Analysis [doi:10.1007/978-3-642-56999-9,doi:10.1007/978-1-4684-6802-1]. The present work quantitatively strengthens this qualitative Turing-completeness: We assign bit-costs to ERC's operational primitives (including partial/multivalued tests) in such a way that any real function incurring polynomial cost becomes Turing-computable in polynomial time, and vice versa. Runtime measurements on implementations in the iRRAM C++ library confirm our theoretical performance predictions.

The Nyström method approximates $A\approx A(\,:\,,I)A(I,I)^{\dagger} A(\,:\,,I)^{\top}=CA(I,I)^{\dagger} C^{\top}$, where $C:=A(:,I)\in\mathbb{R}^{n\times r}$ is a column subset matrix of $A$. When applied to symmetric but indefinite matrices, the Nyström method can fail because the core matrix $A(I,I)$ may severely underestimate the eigenvalues of $A$ and may become (nearly) singular. We address this issue by developing and analyzing an algorithm that carefully chooses $\widehat{M}\in\mathbb{R}^{r\times r}$ in place of $A(I,I)^\dagger$ by solving the two-sided sketched least-squares problem $\min_{M}\|X(A-CMC^{\top})X^{\top}\|_F$, where $X\in\mathbb{R}^{t\times n}$ is a random sketch matrix. We study in detail the cases where $X$ is a Gaussian or a leverage score sampling (LSS) matrix, and show that with oversampling $t>r$ the residual $\|A-C\widehat{M}C^{\top}\|_*$ is comparable to $\min_{M}\|A-CMC^{\top}\|_*$. For the Gaussian sketch, we require $t=\mathcal{O}(r)$ samples; for LSS, we show that $t=\mathcal{O}(r \log r)$ samples suffice for the theoretical guarantee, with the LSS approach carrying the advantage that once a set of $t$ row indices is identified, the approximation requires only $t^{2}$ matrix-entry evaluations to find $\widehat{M}$, given $C$. We illustrate our results with synthetic examples and applications to kernel methods.

We consider the problem of recovering functions defined on the half-line from inexact input information. For the proposed regularizing summation methods based on Laguerre polynomials, we analyze their approximation properties on functions from Wiener classes. We establish conditions under which these methods are stable with respect to small perturbations of the input data and order-optimal not only in accuracy in the weighted $\sup$-norm, but also in the number of Fourier-Laguerre coefficients used.

Let $\Omega\subset\mathbb{R}^d$, $d\ge2$, be a bounded connected domain with boundary of class $C^{1,\alpha}$, where $0<\alpha<1$. For the adjoint Neumann–Poincaré operator $K^*_{\partial\Omega}$, normalised so that its distinguished eigenvalue is $1/2$, let $\lambda_j^+(\Omega)$ denote the upper min–max values on the mean-zero energy space. We prove $\sum_{j=1}^{d}\lambda_j^+(\Omega)\ge \frac{d-2}{2}.$ It follows that $\lambda_1^+(\Omega)\ge \frac{d-2}{2d},$ with equality if and only if $\Omega$ is a ball. In dimension three this proves the $1/6$-conjecture of Miyanishi and Suzuki. The proof uses the coordinate boundary-charge densities induced by uniform applied fields. Their energy Gram matrix is the perfect-conductor polarization tensor $M_\infty$. Positivity of a $2d\times2d$ Gram matrix yields the endpoint Hashin–Shtrikman inequality $|\Omega|\mathrm{Tr}(M_\infty^{-1})\le1,$ and bounds the trace of the compression of $K^*_{\partial\Omega}$ to the applied-field space. Equality in the inverse-trace inequality makes the interior Newtonian potential quadratic, so a converse to Newton's theorem identifies $\Omega$ as an ellipsoid. Equality in the spectral estimate also makes the Hessian of this potential isotropic, which forces the ellipsoid to be a ball.

Selecting the optimal neural-operator prediction during deployment is challenging when high-fidelity reference solutions are unavailable. We demonstrate that under a squared Hilbert-space loss, ranking a finite model library depends strictly on the low-dimensional span of candidate differences, allowing us to score all models simultaneously using a single anchor-based linearized response of the governing equation. This shared physical diagnostic accurately recovered over 99.6% of pairwise preferences and 99.0% of optimal checkpoints across diverse Fourier and convolutional operator libraries for fluid, reaction-diffusion, and wave dynamics. Furthermore, the corrected physical proxy frequently outperformed the best individual candidates, and we establish computable sufficient conditions that rigorously certify exact decisions for strongly monotone discretizations. By exploiting the local dynamical response rather than raw defect magnitude, this framework enables the reliable and highly efficient deployment of scientific surrogates without requiring ground-truth data.

We study a stochastic model for epitaxial thin-film growth driven by spatially rough additive noise in a regime where the noise is singular and regularization via truncation in Fourier space is used to give a meaning to the solution, leading to a vanishing nonlinearity in the limit. In order to study this phenomenon numerically, the nonlinearity is discretized by a spectral Galerkin projection, while time integration is performed with an exponential Euler scheme. For roughness stronger than space-time white noise, we derive strong error estimates in $L^p(\Omega;C([0,T];\mathcal H^1))$ that display explicitly the interaction between the spatial cut-off, the time step, and the decay of the nonlinear current. We also quantify the growth of the truncated stochastic convolution and the corresponding vanishing rate of the nonlinearity. Numerical experiments illustrate the transition from persistent hill formation to noise-dominated dynamics as the roughness parameter increases.

This paper presents a rigorous extension of the explicit, high-order, positivity-preserving, and entropy-stable spectral collocation schemes developed in Upperman 2023 and Yamaleev 2023 for the 3D compressible Navier-Stokes equations to a time-implicit formulation. The time derivative terms are discretized by using the second-order implicit backward difference formula (BDF2) that is well suited for solving time-variable viscous flows at high Reynolds numbers. The nonlinear system of discrete equations resulting from the BDF2 discretization at each physical timestep is solved using a dual time-stepping (DTS) technique. The BDF2 DTS scheme is entropy-stable and positivity-preserving in the pseudotime and provides unconditional stability properties in the physical time. Numerical results demonstrate the efficiency and accuracy of the positivity-preserving BDF2 DTS scheme as compared with its explicit counterpart are presented for supersonic flows with strong shock waves and contact discontinuities.

Sloan iteration raises the convergence order of Galerkin and degenerate-kernel approximations to second-kind integral equations. After quadrature discretization, these become a discrete Galerkin method and a product-integration Nyström method, respectively. How much of this improvement survives quadrature discretization? For zonal integral equations on the sphere, we give a sharp answer by expressing the Sloan error identity in terms of spherical harmonics. In the absence of quadrature, Sloan iteration fully exploits the smoothing of the integral operator. Quadrature can destroy this gain by aliasing unresolved information into low-frequency modes, where further iteration no longer improves the asymptotic rate. This yields a unified tail–aliasing description of these four methods. For positive-weight polynomially exact quadrature and multipliers of exact algebraic order, we derive sharp two-sided worst-case bounds that separate the spectral-tail and aliasing contributions. Their balance yields a depth-dependent recovery–saturation threshold: sufficient overintegration recovers the quadrature-free rate, while below the threshold aliasing determines the sharp order. We also extend the analysis beyond polynomial exactness using Marcinkiewicz–Zygmund stability and Gram-corrected least squares. Numerical experiments illustrate both regimes.

We consider an inverse backscattering problem for time-dependent acoustic waves with a compactly supported penetrable scattering obstacle in unbounded three-dimensional free space. Assuming that dynamic backscattering far field data of scattered waves corresponding to a few time-dependent incident plane waves are available, the goal is to reconstruct the shape of the scattering obstacle. We establish Fréchet differentiability of the time-dependent far field pattern with respect to the shape of the scattering obstacle and give a characterization of the associated temporal domain derivative in terms of its Laplace transform. This characterization is then utilized in an efficient implementation of a regularized Gauß-Newton method for the inverse backscattering problem using convolution quadrature and a boundary element method. Numerical examples demonstrate potentials and limitations of the algorithm.

Tight framelet systems on manifolds provide exact energy preservation and one-pass reconstruction, but their standard semi-discretization relies on polynomial-exact quadrature rules, which are difficult to reconcile with scattered and progressively refined data. We develop nested Marcinkiewicz–Zygmund (MZ) measures on compact Riemannian manifolds as a quantitative alternative. Using dyadic quasi-uniform nested point sets and local partition weights, we establish a sequence of deterministic MZ inequalities for diffusion polynomial spaces. The result has two complementary forms: a fixed-tolerance version, in which the polynomial bandwidth grows dyadically with the level, and a fixed-bandwidth refinement version, in which the MZ tolerance improves as more nested nodes are added. When these measures are used to discretize continuous tight framelets, the MZ tolerance transfers directly to the deviation of the frame bounds from one. Thus quadrature exactness is relaxed with a controlled loss of tightness, and near-tightness improves along the same nested hierarchy. For the fully discrete setting, we develop filter-bank analysis and synthesis procedures, characterize one-pass and canonical reconstruction, and show that the MZ tolerance also controls the conditioning of the frame-operator equation. Fast implementation of filter-bank transforms is also presented. Numerical experiments on the sphere and the flat torus illustrate the resulting multiscale decompositions and reconstruction behavior.

This paper presents a comprehensive, high-order (p-version) Virtual Element Method (VEM) framework for the structural analysis of innovative variable stiffness plates. VEM is particularly suited for complex configurations due to its ability to handle arbitrary polygonal meshes, including curved edges. However, its mathematical formulation may hinder its spread in the engineering community. This work illustrates a formulation with an accessible implementation using well-known FEM notation and integrating at the same time a set of new advanced capabilities. Specifically, both standard stabilized and advanced self-stabilized strategies are adopted. To further improve the robustness of VEM in the presence of variable coefficients, polynomial projections taking into account the coefficients are employed. This approach is referred to as Variable Coefficients-VEM approach (VC-VEM). This unified framework is applied to linear static, free-vibration, and buckling analyses, and validated against analytical solutions and numerical benchmarks. In particular, plates with cutouts and problems featuring high-gradient solutions are investigated, demonstrating that the proposed comprehensive approach provides a flexible and ready-to-implement tool for advanced structural design.

This paper proposes a new technique for computer simulation of linear filters. It allows simulating continuous-time linear filters based on the spectral method for analyzing linear control systems. Applying the spectral method implies that the input and output signals are represented by ordered sets of expansion coefficients in a chosen basis, and the filter itself is specified by a two-dimensional nonstationary transfer function. The described technique is tested on Bessel, Papoulis, and Legendre filters of various orders. For each filter, the corresponding two-dimensional nonstationary transfer function, i.e., the matrix of the linear transformation relating expansion coefficients of the input and output signals, is obtained.

The dimension of a polynomial spline space of bi-degree $(d_1,d_2)$ over a T-mesh $\mathscr{T}$ with the highest order of smoothness $(d_1-1,d_2-1)$ depends on both mesh topology and geometric configurations. Under the assumption that the T-connected components of the T-mesh $\mathscr{T}$ contain no vanishable T $l$-edges, we develop explicit upper and lower bounds of the dimension of the polynomial spline space. By introducing a decoupling technique within the completely non-diagonalizable component (CNDC) of the T-mesh $\mathscr{T}$, we separate tightly coupled multi-vertex constraints and transform global conformality conditions into localized linear equations along each interior large edge. Based on the decoupling technique, a new dimension formula of the polynomial spline space is then presented, and from which sharp upper and lower bounds of the dimension are obtained. The bounds are sharp in the sense that different geometric realizations of T-meshes with the same topology can attain the lower and upper bounds for the dimension of the polynomial spline space. We further prove that the new formula is consistent with Mourrain's homological dimension formula, and a sharper lower bound is obtained by our method.

Shallow networks with prescribed or randomly sampled hidden parameters are widely used as numerical trial spaces, yet their optimal Sobolev approximation power with standard smooth sigmoidal activations in general dimension remains unresolved. We establish the corresponding optimal rates for a class of smooth sigmoidal activations with Schwartz-class derivative decay, including $\tanh$, the logistic sigmoid, and the error function $erf$. We first construct deterministic direction–offset dictionaries with $M$ features such that every $u\in H^k(\Omega)$ can be approximated with error of order $M^{-(k-m)/d}$ in $H^m(\Omega)$ for all $0\le m\le k$. This rate is optimal in the sense of Kolmogorov widths for Sobolev balls. We further prove that dictionaries obtained by independent parameter sampling from any prescribed density bounded away from zero attain the same approximation exponent with high probability, up to logarithmic oversampling. The analysis develops a sigmoidal ridge representation and combines it with deterministic or probabilistic quadrature in direction–offset space while retaining polynomial control of the output coefficients. Numerical experiments across a broad range of dimensions, target regularities, and Sobolev error norms recover the predicted algebraic rates for both deterministic and random feature dictionaries.

Discovering dynamical systems from trajectory data is a central problem in applied mathematics and engineering. Whilst recent advances in machine learning have led to strong progress in data-driven system identification, much less attention has been given to systems with discontinuous dynamics. These systems are nevertheless highly relevant in applications, including climate dynamics and mechanical systems with friction. In this work, we consider the problem of identifying piecewise-smooth dynamical systems directly from trajectory data. Compared with the smooth setting, this requires recovering the governing equations and detecting the switching hyperplanes that separate different dynamical regimes and characterising their behaviour, such as sliding motion. We present a modular framework for discovering such systems by first estimating switching hyperplanes from data and then learning smooth dynamics within each region using geometry-constrained neural networks. The geometry-learning phase is studied from a statistical perspective, analysing the identifiability of the discontinuities and the robustness of the procedure. We also introduce a novel neural network architecture with a prescribed discontinuity set, and provide a theoretical analysis of its approximation properties. The approach is tested on low-dimensional benchmark problems, including dry-friction oscillators and the PP04 climate model for the ice ages.

In this paper, we propose two hierarchical algorithms for charge deposition and electric-field interpolation that apply to both the sparse-grid combination technique (SGCT-PIC) and hierarchical sparse-grid (HSG-PIC) particle-in-cell methods. The two algorithms are inspired by the fast multipole method (FMM) and exploit clusters of particles associated with a directed acyclic graph (DAG) of particle-populated boxes to reduce the number of particle–mesh interactions. The particle–mesh interactions are governed by piecewise-polynomial kernels, so that the associated multipole expansions are exact, requiring neither truncation nor approximation, and are valid in both near- and far-field regions, thereby eliminating the need for multipole-to-local translations. The arithmetic complexity of the charge deposition and field interpolation steps is reduced from $\O(p^d n^{d-1}N)$ to $\O(p^d(N+M))$, where $M=2^{dn}$ denotes the number of full-grid mesh nodes and is typically no larger than the particle population in the considered regime, $M\lesssim N$. Numerical experiments in two-dimensional configurations demonstrate charge-deposition speedups of $8.2\times$–$66.9\times$ for SGCT-PIC and $3.1\times$–$18.8\times$ for HSG-PIC, and field-interpolation speedups of $4.1\times$–$62.6\times$ and $4.2\times$–$13.7\times$, respectively, depending on the particle-per-cell ratio, while preserving the exact particle–mesh interactions. The speedups increase with the particle-per-cell ratio, reflecting the reduced dependence of the hierarchical algorithms on the number of particles and their increasing advantage for large particle populations.

We develop a geometric framework for learning deterministic and stochastic forced Hamiltonian systems with neural networks. Motivated by the Lagrange-d'Alembert principle and the theory of variational integrators, we introduce the notion of a Lagrange-d'Alembert map and establish a $C^r$ convergence theorem for first-order one-step methods. Building on these results, we propose Generalized Forced Hamiltonian Neural Networks (GFHNNs), a class of structure-preserving neural networks obtained by concatenating Lagrange-d'Alembert-Euler maps, and prove a universal approximation theorem for this architecture. We further extend the framework to parameter-dependent systems, leading to Parametric Generalized Forced Hamiltonian Neural Networks (PGFHNNs). By interpreting the multiple Stratonovich integrals appearing in the Stratonovich-Taylor expansion as parameters, the same framework can be applied to stochastic forced Hamiltonian systems whenever information about the underlying Wiener process is available. Our numerical experiments demonstrate that the proposed geometric architectures provide significantly improved long-time stability and accuracy compared to non-geometric residual neural networks, while requiring substantially less training data to achieve a comparable level of performance.

This study demonstrates a rationally enriched Chebyshev (REC) trunk for deep operator network (DeepONet) surrogate models of singularly perturbed and high-Péclet transport problems whose solution profiles are characterized by thin localized boundary or wall layers. The REC trunk combines Chebyshev polynomial dictionary elements with rational dictionary elements constructed using the adaptive Antoulas-Anderson (AAA) algorithm. Over five independent training runs, the resulting REC-trunk DeepONet is evaluated against a vanilla DeepONet and a Chebyshev-trunk DeepONet whose prescribed dictionary consists only of Chebyshev polynomials across three problems whose singular perturbation parameters are diffusion-to-advection ratios: a singularly perturbed scalar boundary-value problem (BVP), the thermal entrance problem with a prescribed wall temperature, and the concentration entrance problem with an absorbing wall. Across the held-out test profiles, the REC-trunk DeepONet improves over the vanilla DeepONet and remains comparable to the Chebyshev-trunk DeepONet in predicting the scalar profile, with its clearest advantage over the Chebyshev-trunk DeepONet appearing when the perturbation parameter lies between $1.00\times10^{-4}$ and $1.78\times10^{-4}$, where it reduces the profile-error metrics by up to $19.5\,\%$ relative to the Chebyshev-trunk DeepONet. In predicting the wall-normal temperature and concentration profiles, the REC-trunk DeepONet reduces the profile-error metrics by up to $60.2\,\%$ and $32.2\,\%$ relative to the vanilla and Chebyshev-trunk DeepONets, respectively, while suppressing artificial near-wall oscillations as the Péclet or mass-transfer Péclet number ranges from $10^{2}$ to $10^{4}$.

In this paper, we develop an evolutionary generalized transferable neural network (Evo-GTransNet) solver for parabolic partial differential equations, formulated as a retained-space fixed-feature Galerkin method of lines. A GTransNet provides the prescribed spatial dictionary, while only the retained output coefficients evolve, thereby avoiding nonlinear training during time integration. To address severe mass-matrix ill-conditioning, we apply a quadrature-weighted truncated singular value decomposition (SVD) to select the numerically resolved trial space, followed by a separate rescaling that makes its basis orthonormal with respect to the assembly-quadrature mass inner product. Rank truncation modifies the approximation space, whereas the subsequent orthonormalization changes only its coordinate representation and preserves the retained discrete functions in exact arithmetic. The resulting semidiscrete coefficient system has an identity mass matrix, and we establish a semidiscrete energy law for symmetric linear parabolic problems. With the implicit midpoint scheme for time discretization, we further prove the contractivity of the method and derive a conditional fully discrete error estimate in which the error is controlled by the retained-space approximation error and the consistency defects. Numerical experiments demonstrate second-order temporal convergence using repeated feature samples together with separate assembly and validation quadratures, and quantify the accuracy of the retained space in the presence of severe raw-mass ill-conditioning. For the high-frequency and multiscale benchmark problems considered here, GTransNet achieves the smallest mean validation errors among the tested fixed-feature dictionaries at the same nominal output dimension.

We present a macroscopically conservative, rank-adaptive method for solving a hybrid Vlasov-Fokker-Planck (VFP) equation in which the ions are treated kinetically, and the electrons are treated as a fluid. Solving this system poses several coupled computational challenges: the curse of dimensionality, conservation of macroscopic quantities, stiffness arising from the multiscale nature of the model, and structure preservation. To address these challenges, we combine several established methods, each targeting a specific difficulty, into a unified framework for solving this nonlinear system. Specifically, we employ the recent Reduced Augmentation Implicit Low-rank (RAIL) method for rank adaptivity in velocity space to combat the curse of dimensionality, high-order implicit-explicit time stepping to handle the multiscale stiffness, and the Local Macroscopic Conservative (LoMaC) procedure for conservative truncation of the solution. The RAIL and LoMaC methods are extended from Cartesian to cylindrical coordinates. The resulting framework accommodates implicit rank-adaptive time integration while remaining conservative and leveraging structure-preserving discretizations. We verify the scheme on a suite of test problems and simulate a standing shock to demonstrate the importance of conserving the macroscopic quantities.

For given $m \times n$ data matrices $X, Y$, we investigate the Normal Procrustes Problem—the least squares optimization problem that aims to minimize $\|AX-Y\|_F^2$, where $A$ is constrained to be a normal $m \times m$ matrix. As far as the author of this article is aware, no other method that attempts to solve the Normal Procrustes Problem exists in the literature; we thus propose what is, to our knowledge, the first such method. We, furthermore, adapt our approach to address the Real Normal Procrustes Problem, where $A$ must be real. In our treatment of these problems, we first reduce our complex and real objective functions to be purely optimizable over the Riemannian manifolds of the unitary and real orthogonal matrices, respectively. This reduction enables us to apply techniques in Riemannian manifold optimization to approximate solutions to both. The Closest Normal Matrix and Real Closest Normal Matrix Problems are both special cases of their respective Procrustes Problems (by setting $m = n$ and $X = I_m$) and have been previously studied in the literature. Our approach thus recovers a novel Riemannian optimization method for approximating solutions to both these problems. We further numerically test the performance of our method across all such problems (including against previously developed algorithms on the Closest Normal Matrix Problems) and obtain competitive residuals and favorable scaling in wall-clock time.

We implement Ozaki Scheme II (residue number system + Chinese remainder theorem), which reduces multiple-precision dense matrix multiplication to a sequence of low-precision, high-throughput integer or floating-point GEMMs on CPUs. Two backends are built on top of a shared CRT reconstruction stage: (a) exact INT8 x INT8 -> INT32 tile products on Intel AMX, and (b) binary64 DGEMM, the CPU construction of the original Ozaki Scheme II paper. On a two-socket Xeon Gold 6526Y (Emerald Rapids, 32 cores), we evaluate significand precisions of 53-2048 bits and matrix dimensions N = 256-8192. The results are always within 1 ulp of a high-precision MPFR reference (essentially correctly rounded), while running up to 167x faster than a naive MPFR matrix product, up to 588x faster than BNCmatmul's Strassen multiplication, and 9-78x faster than Ozaki Scheme I (FP64 slicing + OpenBLAS DGEMM). The break-even point between the two backends is approximately N = 2048: below it, the binary64 backend wins thanks to its smaller number of moduli; above it, the AMX-INT8 backend wins as the GEMMs dominate. We further port the implementation to AArch64 (NVIDIA GB10: Cortex-X925 x 10 + Cortex-A725 x 10). Since this machine lacks SME/SME2, the INT8 kernel uses the SMMLA matrix-product instruction of the SVE2 i8mm extension. We obtain an exact INT8 GEMM sustaining 6.5 TOPS and, still within 1 ulp across all conditions, speedups of 14-89x over BNCmatmul's Ozaki Scheme I (OpenBLAS-linked routine) and 6-19x over a fairness-adjusted OzI-best variant. The paper also includes a tutorial introduction to Ozaki Scheme II (Section "Introduction to Ozaki Scheme II") and a quantitative explanation of why this seemingly GPU-specific technique is fast on CPUs as well.

Quadratic manifold-based model order reduction offers a viable pathway to circumvent the limitations of linear subspaces for linear control systems characterized by slow Kolmogorov $n$-width decay. However, a system-theoretic framework for constructing such quadratic approximations remains absent from the literature. This paper presents a system-agnostic, optimization-free framework for the direct construction of quadratic projection matrices. We prove that the synthesized reduced-order model matches the nonlinear moments of the full-order system and preserves its exact center manifold mapping, thereby ensuring asymptotic tracking of steady-state outputs under specific input classes. Numerical results on transport-dominated benchmark problems, namely, the one-dimensional damped wave and advection equations, show that the proposed framework achieves high-fidelity trajectory reconstruction within a significantly reduced-dimensional state space, yielding substantial online computational savings.

We consider the task of simulating the real-time dynamics of weakly interacting fermionic systems. In particular, we focus on computing the expectation value of a local observable $A$ at time $t$. By analyzing the convergence of the perturbative expansion in the interaction strength $\lambda$ for the Heisenberg-picture observable, we propose a polynomial-time algorithm for estimating this expectation value in the weakly interacting regime $\lambda |t|^{2D+1}=\mathcal{O}(1)$, when the Hamiltonian is geometrically local on a $D$-dimensional lattice. Importantly, this condition is independent of the system size. If the goal is instead to approximate the time-evolved observable in normalized Frobenius norm, we extend the convergence regime to $\lambda |t|=\mathcal{O}(1)$ with quasi-polynomial runtime. When the non-interacting part exhibits Anderson localization, our polynomial-time algorithm can be extended up to $\lambda |t|=\mathcal{O}(1)$, modulo polylogarithmic factors. Our algorithm brings together ideas from continuous-time QMC, diagrammatic QMC, and Majorana Propagation, but with a new Heisenberg-picture operator-growth analysis that makes the sampling complexity rigorously controllable. This leads to provably efficient classical algorithms in regimes where the interaction is weak enough that the sampling variance remains bounded independently of system size. Together, these results identify broad regimes in which weak interactions, locality, and localization can be leveraged to make real-time fermionic dynamics classically tractable.

In this study, a two-dimensional second-order macroscopic traffic flow model is proposed to capture the complex dynamics of disordered traffic by explicitly incorporating coupled longitudinal and lateral interactions. The framework extends the classical continuity equation to two spatial dimensions and introduces acceleration equations consisting of self-driven, interaction-induced, and road-boundary-related components. The macroscopic longitudinal acceleration is formulated based on the Full Velocity Difference Model (FVDM), while lateral dynamics are governed by interaction principles consistent with the Optimal Velocity Model (OVM). The explicit inclusion of road-boundary effects enables the representation of vehicle confinement and space-sharing behavior that are essential in lane-free and weakly lane-disciplined traffic systems. The model is examined through a series of numerical experiments in which longitudinal and lateral dynamics are analysed individually as well as simultaneously within a coupled two-dimensional setting under a range of initial conditions, including transitions between free-flow, medium congestion, and heavy congestion, and different lateral density configurations. The simulations demonstrate the model's ability to reproduce key traffic features such as shockwave propagation, lateral dispersion, vehicle rearrangement, and stable density evolution across the road width. The numerical scheme, based on upwind and Lax-Friedrichs discretisations, ensures stable solutions of the coupled partial differential equations. Overall, the proposed framework provides a robust macroscopic description of disordered traffic and offers a consistent basis for analysing two-dimensional vehicular flow dynamics.

Designing compact, high-efficiency thermal architectures requires resolving the competing demands of solid conduction, fluid convection, and flow resistance within highly constrained physical envelopes. Density-based topology optimization provides a systematic framework for synthesizing these coupled layouts, yet the reproducibility and numerical stability of the resulting designs depend critically on the underlying discrete solvers and adjoint sensitivity mechanics. In this work, we present a transparent, self-contained two-dimensional finite-volume formulation on a staggered Marker-and-Cell grid for conjugate heat transfer governed by design-dependent energy transport coupled to Stokes–Brinkman or Darcy flow at fixed solid volume. To prevent spurious artificial thermal sources in porous, weakly compressible Brinkman domains, the discrete advection operator is constructed to satisfy the identity $\mathbf{u}\cdot\nabla T=\nabla\cdot(\mathbf{u}T)-T(\nabla\cdot\mathbf{u})$ cellwise, ensuring that uniform temperature fields remain exact discrete nullspaces even under inexact continuity satisfaction. Reverse-mode derivatives are evaluated via the implicit function theorem rather than unrolled iterative loops, yielding exact discrete adjoints with bounded memory requirements. The discrete operators are systematically validated through the method of manufactured solutions and directional Taylor remainder tests. Four representative thermofluid design benchmarks are optimized using a projected-gradient scheme with $\beta$-continuation, wherein candidate iterates are accepted and published only upon satisfying rigorous, predeclared gates on residual convergence, mass conservation, volume feasibility, and numerical finiteness. The resulting formulation provides an inspectable, deterministic reference stack for verifiable conjugate thermofluidic topology optimization.

We present a high-order radial basis function-generated finite difference (RBF-FD) method for partial differential equations on moving manifolds $\mathcal M(t)\subset\mathbb R^3$ of co-dimension one. Our method builds on the tangent-plane formulation of surface RBF-FD and combines Lagrangian and Eulerian treatments: the manifold and material derivative are evolved in a Lagrangian fashion, while the remaining surface differential operators are reconstructed on the instantaneous point cloud and stabilized, when necessary, by a quasi-analytical hyperviscosity formulation. Lagrangian marker drift is handled by adaptive rearrangement using a compact global parametric model of the moving surface, which also supplies accurate normals and geometry-based quadrature. After rearrangement, we reconstruct the multistep history by backward semi-Lagrangian tracing and interpolation before resuming the Lagrangian time discretization. We exploit temporal coherence through three update strategies: defect correction for the local RBF-FD weights, a curvature-based update of the hyperviscosity coefficients between spectral recomputations, and a global defect-correction iteration that reuses an incomplete LU (ILU) factorization before preconditioned generalized minimal residual (GMRES) iterations. Finally, we enforce the prescribed global mass balance through a scalar projection based on the evolving surface quadrature. Numerical experiments demonstrate high-order convergence, conservation to roundoff in source-free problems, stable long-time integration, and substantial savings from the proposed update strategies.

This paper presents a high-order rank-adaptive implicit integrator for the tensor solution of high-dimensional diffusion equations. We extend the 3D version of this method from the Tucker decomposition to higher dimensions using the hierarchical Tucker (HT) decomposition, since the storage complexity for the Tucker decomposition increases exponentially with the number of dimensions $d>3$. The HT format avoids this issue by decomposing the solution according to a binary tree consisting of bases for each dimension and core tensors which connect the bases. Spectral methods are considered for spatial discretization, and diagonally implicit Runge-Kutta methods are considered for time discretization. At each stage of the Runge-Kutta method, the bases computed at the previous stages are augmented to predict the upcoming basis and construct projection subspaces. By projecting onto these enriched subspaces, the bases and cores can be updated in a sequential manner going up the tree from leaf-to-root. Unlike the 3D Tucker method which has a single core tensor, the HT method also updates the intermediate core tensors. Numerical experiments demonstrate that the method observes high-order accuracy, and test how well the integrator captures the solution rank for various sets of time-dependent diffusion coefficients.