Papers

Recently, efficient sweep-based low-rank iterative solvers have been developed to reduce the computational costs of steady-state solves and implicit time integration of the radiative transfer equation (RTE). In particular, angular sampling based on the discrete empirical interpolation method (DEIM) or on error-indicated randomized greedy selection enables efficient, non-intrusive reuse of existing transport sweep implementations. However, randomized-greedy low-rank source iteration with diffusion synthetic acceleration (SI–DSA) enforces uniformly tight accuracy tolerances across all source iterations, thereby oversolving the early iterations, and it does not explicitly reuse angular bases from preceding iterations. More critically, like its full-rank counterpart, this approach may lose effectiveness in multiscale problems with strong material discontinuities. We address these limitations through two complementary developments. First, we introduce adaptive accuracy control into randomized-greedy low-rank SI–DSA, reducing intermediate ranks and computational cost while maintaining the accuracy required for convergence. A two-stage hybrid sampling strategy can provide further acceleration by switching to DEIM-based sampling near convergence, once the preceding angular basis becomes sufficiently representative. Second, we develop an inexact flexible GMRES–DSA solver for the integral scalar-flux reformulation, in which sweep-based low-rank approximations are confined to the construction of the right-hand side and the matrix–vector products. The scalar-flux Krylov basis is not compressed, thereby preserving Arnoldi orthogonality in exact arithmetic and allowing established inexact-GMRES analysis to guide adaptive accuracy control of low-rank operations. We present numerical results on benchmark problems that demonstrate the effectiveness of the proposed methods.

Inverse inequalities are essential for controlling derivative growth and ensuring stable sampling in kernel approximation. Existing whole-space and single-scale estimates provide an established analytical basis, but their extension to multiscale analytic spaces on finite intervals requires uniform control of boundary effects and cross-scale cancellation. This paper establishes unweighted interval \(L^2\) inverse inequalities for spaces generated by primitives of kernels with positive exponential spectra. By retaining exterior diagonal energy and controlling off-diagonal interactions through one-sided estimates with endpoint-distance weights, we derive uniform coefficient bounds for normalized high derivatives. A synthesis-operator perturbation argument extends these bounds to positive spectral mixtures, while adjacent-derivative comparison and interval interpolation yield an inverse bound of order \(h_*^{-1}\), where \(h_*\) is the smallest scale. Under geometric scale separation and within-scale centre separation, the constant is independent of the number of scales and centres and their distances from the endpoints, allowing endpoint centres and coincident centres across scales. The coefficient lower bound approaches the optimal limit \(1/2\). Two-centre constructions establish the sharpness of the scale exponent for every fixed kernel in the class and, for a single exponential spectrum, a sharp deficit of order \(q^{-1/2}\) from \(1/2\) for the best uniform coefficient lower bound, where \(q\) is the derivative order. Numerical illustrations for Cauchy and logistic kernels demonstrate coefficient stability in several multiscale configurations. These results extend interval inverse estimates to multiscale analytic kernel spaces with endpoint centres and provide quantitative conditions for deterministic sampling stability.

We develop a mesh-free neural JKO scheme for advection-reaction-diffusion equations with a gradient-flow structure in the Hellinger-Kantorovich (HK) geometry of unbalanced optimal transport. Each update is parametrized by a spatial map and a mass-changing factor, allowing spatial redistribution and local mass creation or loss to be treated jointly within a single variational step. Their cone action bounds the squared HK distance from above, yielding a sufficient condition for discrete energy dissipation through comparison with the identity pair. Minimizing the pair objective over all admissible pairs recovers the exact JKO minimum when the source and a minimizer have positive densities. We establish existence and mass bounds for JKO minimizers and, under additional assumptions, obtain positivity and regularity together with a discrete Euler-Lagrange equation and a metric-dissipation identity. The self-consistent chemical potential is then nonincreasing along an optimal map. There exist parametric pairs whose endpoint densities and objective values converge to those of an exact JKO minimizer, provided a regular-pair approximation hypothesis holds. Finally, we show that a primal-dual gap controls objective suboptimality and, for Boltzmann entropy, the $L^1$ density error, assuming exact-step regularity, positive-semidefinite interactions, and global dual feasibility. Numerical experiments examine pointwise agreement with the PDE, energy dissipation, and the roles of transport, reaction, and fully implicit interactions.

When a feasible plan is obtained whose quadratic transport cost is known to be close to the optimal cost, we try to determine how close the feasible plan is to the true optimal map under mild density and moment conditions. The source and target may be unbounded or have non-compact supports. We show that for a source with density in $L^p$ and an $n$-th moment, let $s=1-1/p$ and $\eta=sn/[n+(d-1)s]$. Under the target constraint $\int |y|^m[\log(e+|y|)]^\beta\,d\nu\leqslant1$, the worst-case class-uniform conditioning has a sharp modulus comparable to \[ t^{\eta(m-2)/[m(1+\eta)-\eta]} [\log(e/t)]^{-\beta(2+\eta)/[m(1+\eta)-\eta]}, \] whenever $t>0$ is small. This holds for $m\geqslant2$ and $\beta\geqslant0$. When $m=2$, every $\beta>0$ gives a sharp logarithmic modulus, whereas the $m=2,\,\beta=0$ extreme has no vanishing modulus. The matching lower bound is verified for both the squared map error and barycentric projection error. The sharp map error bounds allow us to directly derive quantitative stabilities of common source Brenier maps without passing through potential estimates, thereby avoiding loss of information. These stabilities strictly improve the corresponding results of Delalande-Merigot \cite{delalande2023quantitative} and Letrouit-Mérigot \cite{letrouit2026gluing} under identical or weaker settings, particularly no convexity of the source support or bounds of the source density are assumed.

We develop and analyze a B-spline scaling-function Galerkin method for the tempered fractional partial differential equation governing European and American option prices under the CGMY process. The spatial generator couples exponentially tempered left- and right-sided Riemann–Liouville derivatives of order $Y\in(0,2)$, whose nonlocality and asymmetry present discretization challenges. We derive a corrected Fourier-symbol construction of the tempered connection-coefficient matrices whose interior blocks are Toeplitz, maintain exact fractional scaling, and enable $\mathcal{O}(N\log N)$ matrix–vector products via the fast Fourier transform. In the regime $Y\in(1,2)$, we prove full-line and finite-domain Gårding inequalities, coercive whenever the risk-free rate is positive, via zero-extension to the full-line symbol. On this foundation, we establish elliptic-projection, semidiscrete, and fully discrete error estimates with spatial rate $\mathcal{O}(h^{p-Y/2})$, alongside an unconditional shifted-energy stability estimate for the forced Crank–Nicolson scheme with mesh-independent constants. For positive rates the admissible shift vanishes, governing the plain scheme exactly as implemented. A single backward solve prices an entire strike strip at fixed maturity with smooth, analytically differentiated Greeks. The same discretization prices American puts through a projected linear complementarity solve, with early-exercise premiums and exercise boundaries validated against a resolved deterministic reference and a Longstaff–Schwartz Monte Carlo benchmark. Calibration to S&P~500 option surfaces across low-volatility, crisis, and recovery regimes yields implied-volatility errors near two vol-points on well-identified surfaces, while an accuracy-matched comparison with a stabilized Grünwald–Letnikov baseline delineates precisely when the high-order variational structure pays off.

Approximating linear operators on high-dimensional spaces using only forward operator-vector product queries is a common but challenging task in scientific computing and machine learning arising when access to the adjoint action is experimentally unavailable or difficult to construct computationally. Existing query-efficient methods with polylogarithmic query complexity in the error tolerance rely on self-adjointness or access to the adjoint operator, while existing adjoint-free methods based on global low-order regularity have polynomial query complexity. We show that a priori local row-space information can be exploited for query-efficient, adjoint-free operator approximation, and that this information can be supplied by local high-order regularity estimates. Our approach first extends the fixed-sparsity matrix approximation algorithm of (Amsel et al., SIMAX, 2026) to approximating a matrix by one whose rows lie in prescribed subspaces, providing a general template for exploiting local row-space structure from forward queries alone. We then apply this framework to Hilbert-Schmidt integral operators with asymptotically smooth kernels. We establish asymptotic smoothness for Green's functions of canonical second-order elliptic PDEs in $d\leq 3$ dimensions including Poisson's equation in balls and constant coefficient equations on tori. In $d$ spatial dimensions, our algorithm uses $\mathcal{O}\big(\log(\varepsilon^{-1})\log^{d}(\varepsilon^{-1}\delta^{-1/2})\big)$ non-adaptive randomized queries to produce an $\varepsilon$-accurate approximation in Hilbert-Schmidt norm with probability $1-\delta$. This implies the error decreases as $\mathcal{O}\big( \exp(- c q^{1/(d+1)}) \big)$ in the query count $q$, which we observe in numerical experiments using tens to hundreds of queries to approximate solution operators of one-dimensional variable-coefficient second-order elliptic equations.

In this paper, we propose an alternating numerical scheme for constructing quasi-periodic solutions to a class of nonlinear partial differential equations (PDEs) with bounded perturbations, including the one-dimensional nonlinear Schrödinger (NLS) equations and nonlinear wave (NLW) equations under periodic and Dirichlet boundary conditions. A key observation concerning the resonant $Q$-equations allows us to establish directly a local diffeomorphism between the amplitudes and the drifted frequencies. This renders the proof considerably cleaner, since it avoids the Birkhoff normal form and the associated coordinate transformations required in classical KAM theory for PDEs.Unlike in the finite-dimensional setting, the lattice boxes must also include the spatial index, so that the multiscale induction has to be carried out in one additional dimension. This introduces two difficulties. First, the measure estimate for the excluded set of “bad” frequencies must account for the spatial index. For both the NLS and NLW equations, the linear frequencies and their differences are integers or close to integers; this allows the spatial dependence to be controlled through the modulus of the temporal index, so that the excluded measure remains under control. Second, a singular set may contain sites with the same temporal index but widely separated spatial indices, which obstructs the “inversion implies localization” process. For the one-dimensional NLS and NLW equations, however, such sites occur only in pairs with opposite spatial indices. They can therefore be separated by small boxes, after which the global inverse is obtained from the local inverses via the resolvent identity. Numerical experiments illustrate the scheme for nonlinear traveling and standing waves, exhibiting both the computed quasi-periodic solutions and the convergence of the iteration.

Learning constitutive laws of complex fluids from velocity data (i.e. indirect observations) is a PDE-constrained inverse problem in which expressive neural parameterizations risk breaking the well-posedness of the forward physics model. We address this tension by learning, rather than the constitutive law itself, the dissipation potential: a scalar function whose convexity, frame-indifference, and dissipativity propagate into the underlying continuum mechanics and guarantee the structural properties needed for a well-posed and generalizable forward PDE. We introduce ICNNE, an input-convex neural architecture that exactly enforces convexity, frame-indifference, and evenness in the second strain-rate invariant via symmetrization, alleviating the numerical instabilities that occur at small strain rates in standard formulations. By weakly enforcing zero gradients in the origin of the potential via penalization, ICNNE also satisfies the dissipativity property. The loss objective and its gradient are computed using a combination of finite element and neural network methods, coupling the Firedrake and PyTorch libraries. We evaluate on four problems spanning compressible and incompressible regimes: compressible Navier-Stokes, Herschel-Bulkley yield-stress flow, Hibler's viscous-plastic sea-ice model, and discrete-element-method data with no known constitutive law. We show that the learned potentials (i) recover ground-truth physics where it is known, (ii) transfer accurately to geometries unseen during training, and (iii) succeed in regimes where unstructured methods diverge. These results suggest that embedding mathematical structure into the parameterization, rather than into the loss, is a robust path to learning physics from indirect data.

Systems of equations arise in a wide range of scientific and engineering applications. The present work focuses on solvers for general systems of equations, including but not limited to those arising from partial differential equations. These systems can be broadly categorized into linear and nonlinear problems. For large linear systems, iterative solvers are generally preferred over direct methods due to the latter's superlinear growth of computational costs. Although convergence theory is well-developed under certain assumptions on the coefficient matrix, many classes of systems still pose open challenges. These difficulties become even more severe for systems of nonlinear equations, where nonlinear solvers typically rely on repeated linearization. For example, Newton's method may even converge quadratically near the solution; it can also converge slowly or diverge when the initial guess is not chosen appropriately. A wide range of solvers with diverse variants and hyperparameter settings exists, and the development of efficient and robust iterative methods remains an active area of research. Recently, machine learning (ML) techniques have been applied to enhance the efficiency of classical iterative methods while preserving their interpretability and reliability. We refer to these ML-enhanced iterative methods as hybrid iterative methods, in the sense that they combine classical iterative methods with ML. This paper provides a comprehensive overview of state-of-the-art approaches to constructing hybrid iterative methods for systems of both linear and nonlinear equations, while also discussing open challenges and outlining potential directions for future research.

We propose the safeguarded Anderson-accelerated reflection (SAAR) method, a two-level iterative method for the continuous Sylvester matrix equation $AX+XB=C$. The core of the method is a fixed-point reformulation: at each outer iteration, the associated linear matrix equation is decomposed into independent column-wise linear systems, which are solved approximately by a reflection iteration. For nonsingular \(A\) and \(n\ge 2\), we prove the convergence of the reflection inner iteration by a spectral-radius argument and extend the result to the matrix equation setting. Since only finitely many inner iterations are used in practice, the overall scheme is analyzed as an inexact fixed-point iteration, and a sufficient condition for linear convergence of the unaccelerated inexact iteration is established in the Frobenius norm. On top of this scheme, SAAR incorporates Anderson acceleration into the outer iteration, while a residual-based safeguard rejects Anderson candidates whose Sylvester residual is larger than that of the standard inexact reflection iterate, thereby enhancing numerical stability. Numerical experiments indicate that SAAR substantially reduces the number of outer iterations and often the CPU time, and it improves empirical robustness with respect to the inner tolerance for the tested problems.

The development of efficient iterative schemes for approximating the inverse and generalized inverse matrices is important in my research fields. To this end, a family of fourth-order matrix iterative schemes for approximating the inverse and generalized inverse matrices is established by incorporating an appropriate weight function into the scalar Heun method. Moreover, a fifth-order matrix iterative scheme is obtained as a particular member of the proposed family. Under suitable conditions, the theoretical order of convergence is determined for the proposed family, whose stable members are identified through dynamical analysis. To assess their performance and verify the consistency of the theoretical results, some numerical tests are performed on matrices of different sizes.

Gaussian elimination is the oldest and most popular method for solving a linear system. Its numerical stability for a given matrix is controlled by the growth factor, a measure of how large entries can become during elimination. In this work, we prove a number of new results regarding the growth factor of random matrices. First, we provide tight estimates for the growth factor of a matrix preconditioned by a Haar orthogonal matrix without pivoting and characterize the asymptotic distribution of the growth factor of Gaussian matrices without pivoting. Second, and most notably, we prove that the growth factor of an $n \times n$ Gaussian matrix under partial pivoting is rarely much larger than $\sqrt{n}$, resolving an old conjecture of Nick Trefethen. The same techniques used to prove this conjecture also provide an improved smoothed analysis of the growth factor under partial pivoting.

Neural solvers for partial differential equations (PDEs) can approach an accurate solution while their parameters grow without bound. In such cases, the limiting solution may have no finite representation in the chosen model, leaving the best loss unattained. Our analysis connects missing limits in deep neural tanh- networks to unbounded hidden parameters or increasingly redundant neurons. For a class of models built from translated kernels, we describe the missing functions and recover them by adding kernel derivatives to the model. This completion makes the best approximation attainable under standard assumptions. Numerical studies follow the associated parameter growth and explore how completion affects PDE optimization.

This paper introduces a new extrapolation framework within the Calculus of Variations. We develop a robust quantitative theory of the $\Gamma$-convergence approach of De Giorgi by combining it with Gagliardo's work on relative completions in interpolation theory. A central new tool is the extrapolation of \emph{variational $K$-functionals}, which associates to an arbitrary family of functionals a dynamic penalty method that simultaneously characterizes their $\Gamma$-limits and yields explicit rates of convergence for minimizers. As a technical foundation we prove compactness extrapolation theorems also involving Gagliardo completions. To illustrate the mechanism, we apply the technique to variational models arising in image processing, notably the Rudin–Osher–Fatemi (ROF) and Chan–Esedoglu models. Our main result in this setting is a contribution to the inverse regularity problem associated to discrete approximations of the ROF model: we show that the required regularity of the observed image that guarantees a prescribed optimal rate of convergence is measured in terms of the decay of variational $K$-functionals. Surprisingly, our characterization is explicitly captured by using the recently introduced spaces of Brezis, Van Schaftingen, and Yung. Moreover, our methods unify, and considerably extend, classical pointwise convergence and $\Gamma$-versions of the Bourgain–Brezis–Mironescu and Maz'ya–Shaposhnikova formulae, as well as related results in a variety of settings that involve BMO-type functionals, heat kernels, Sobolev spaces on metric measure spaces, ... By bridging De Giorgi's variational convergence with Gagliardo's completion through variational $K$-functionals and extrapolation of compactness, this work establishes a quantitative approach to the Calculus of Variations.

In this work, we study transparent boundary conditions for the time-fractional heat equation on a metric star graph. This equation is the standard continuum model for subdiffusive heat and mass transport in disordered and porous branched media. Using the Laplace transform to solve the exterior branch problems, we demonstrate that the non-local vertex operator observed by the incoming bond is an Abel-type convolution kernel of order $\beta/2$, which generalizes the classical half-derivative kernel recovered at $\beta=1$. Our main result is that the diffusion-coefficient sum rule $\kappa_1=\kappa_2+\kappa_3$, known to eliminate thermal backflow for the classical heat equation, remains unchanged for every fractional order $\beta$: the transparency of the junction depends only on the network's diffusivities and not on the memory of the transport law. We rigorously prove this and uniqueness for the coupled vertex problem. We confirm our results numerically at the level of the underlying coupled partial differential equation system with an implicit L1 finite-difference scheme. We compare the star-graph solution to a reference solution on an unbranched line with matching diffusivity.

We study the surface Stokes-Cahn-Hilliard system governing tangential flow and phase separation in viscous fluid membranes. The equations are posed on a closed, stationary $C^2$ hypersurface $\Gamma\subset\mathbb{R}^{d+1}$, $d=2,3$. For constant viscosity and mobility and a class of regular double-well potentials, we prove global existence and uniqueness of weak solutions for initial phase fields in $H^1(\Gamma)$, with the velocity chosen orthogonal to the Killing fields. These solutions satisfy an energy equality and conserve mass. We establish Lipschitz continuous dependence in energy norms on every finite time interval, including for initial data with different means, and derive higher regularity results for $d=2$ under additional smoothness assumptions on the surface, potential, and initial data. We also develop a fully discrete unfitted finite element method that conserves the total phase and satisfies an exact discrete energy identity. The method combines a penalty on approximate Killing fields with projections in the transport and capillary terms to control these modes while preserving the energy balance. Numerical experiments examine the identification of Killing fields and illustrate phase separation on a torus and a liposome with experimentally informed material parameters.

Let all pairwise distances of an arbitrary configuration of points in $\mathbb{R}^{\ell}$ be known exactly. We prove that every second-order stationary point of the complete-graph squared-stress objective over configurations in $\mathbb{R}^{k}$ has zero residual whenever $k > \ell$. The result holds for every number of points, without a genericity assumption, and allows repeated points and degenerate ground-truth configurations. It resolves the one-extra-dimension conjecture for complete, uniformly weighted data. The proof combines positivity of the residual stress, polar normalization of the candidate configuration, a row-wise rank-nullity argument, and propagation of zero curvature through the full Hessian. The theorem, including the passage from derivatives of the original quartic objective to its matrix formulation, has been formalized in Lean 4. The accompanying verification records include a kernel check with zero trust level and a transitive axiom audit.

We present an algorithm for handling the small cell problem in Cartesian cut cell grids. The advantage of this variation is that it preserves the base scheme's steady state exactly. We show several computational examples using finite volume schemes in two space dimensions.

The performance of stationary iterative methods is memory-bound, due to the low arithmetic intensity. So, we propose a mixed-precision technique to reduce the amount of data movement in matrix-splitting iterative methods. We demonstrate both theoretically and experimentally that this technique will generally converge at a similar rate as an exact-precision iteration, depending on the precision used and the conditioning of the $M$ matrix. Furthermore, we demonstrate that iterative refinement can be used to refine the solution to full double-precision accuracy. Using this approach, we mixed half and double precision to achieve an average speedup over uniform double precision of $1.60\times$ for scalar and block Jacobi and $1.25\times$ for scalar and block Gauss-Seidel. Furthermore, we reduced the precision in the SSOR iteration used by OVERFLOW, a NASA code for compressible fluid dynamics simulations, and achieved a $1.14\times$ speedup to the overall application.

Isogeometric analysis (IGA) solves partial differential equations accurately on exact NURBS geometry, whereas neural solvers are mesh-free but often orders of magnitude less accurate and typically trained by non-convex optimization without error control. We propose IGA-KAN, which uses local Kolmogorov-Arnold networks, fitted in closed form, to improve the IGA solution instead of replacing it. An IGA Galerkin solve produces u_h; on every knot-vertex patch a Kolmogorov-Arnold ridge model is fitted to the strong form of the equation, the exact boundary data and u_h, and the models are blended by IGA hat functions. With fixed inner functions the fit is one batched linear least-squares problem, without optimizer, learning rate or initialization. An a posteriori safeguard, motivated by a maximum-principle bound, decides where local models are used, keeping the IGA solution elsewhere. On eight benchmarks with exact solutions, five from the literature and one also posed on a domain fitted to a brain slice from MRI, the method reduces the error of IGA, at an unchanged number of Galerkin unknowns, by factors of 4.2 to 90 in L^2 and 4.1 to 220 in H^1 on the reference meshes, and its L^2 error is 6 to 6x10^4 times smaller than that of the best Kolmogorov-Arnold network trained from scratch on the same equations with a fixed budget. In an inverse problem it recovers an unknown constant source from one noise-free observation 167 times more accurately than IGA. The gain is attributed to the superconvergence of local averages of the Galerkin solution.

We study deterministic and randomized midpoint discretizations of Langevin dynamics for a target $\pi \propto e^{-V}$, where $0 \prec \alpha I\preceq\nabla^2V\preceq\beta I$ and $\kappa=\beta/\alpha$. To achieve $\sqrt\alpha\,W_2\leqslant\varepsilon$, we show that deterministic Heun uses at most $\widetilde O(\kappa^{4/3}d^{1/3}\varepsilon^{-2/3})$ gradient queries, and underdamped exponential midpoint uses $\widetilde O(\kappa^{5/4}d^{1/4}\varepsilon^{-1/2})$. The proofs exploit cancellation at stationarity and smoothing using techniques from Malliavin calculus, outperforming previous upper bounds based on standard couplings. At bounded condition number, a lower bound matches the $d$ and $\varepsilon$ powers of both deterministic methods. To contrast, for the randomized midpoint methods and Poisson midpoint with at least two grid points (both overdamped and underdamped variants), a simple Gaussian calculation yields a lower bound $d^{1/3}\varepsilon^{-1/3}$ to get an $\varepsilon$-close sample despite starting at a benign initialization. This shows surprisingly that in high dimensions, deterministic discretizations can outperform their random counterparts.

Physics-informed machine learning as recently emerged as a new paradigm for the resolution of partial differential equations, including functionalities for parametric modeling and inverse design, using neural networks to represent the solution fields. However, several difficulties have been reported concerning the training of such networks, resulting in slow convergence, low accurate solutions and sensitivity to the choice of hyper-parameters. In order to better understand the mechanisms underlying these difficulties, the present work proposes a rigourous comparison of neural networks and B-Splines representations within the same physics-informed learning framework. Four problems with closed-form solutions are numerically studied to establish a rigorous analysis of the accuracy and efficiency of the two representations, in terms of asymptotic convergence rates, sampling requirements, loss function minimization and arithmetic precision needed.

In this work, we propose and analyze low-order nonconforming face-based finite element methods for elliptic problems on complex hybrid three-dimensional meshes, composed of arbitrary polytopal cells. The focus is on scalar problems, but the approach can be extended in a tensor manner to design generalizations of the classical Crouzeix–Raviart and Rannacher–Turek methods for mixed formulations. We adopt a nonparametric construction of the local approximation spaces. We demonstrate that it avoids the degradation of approximation properties commonly induced by parametric mappings on distorted meshes. Using these approximation spaces, we define finite element schemes for very general cells. We further demonstrate that, on strongly degraded meshes, an approach inspired by discontinuous Galerkin methods enables removing consistency errors, even in the presence of potentially curved faces. Explicit examples of the construction of the polynomial space are given for prismatic (of triangular basis) and pyramidal (of quadrilateral basis) cells. The resulting schemes are shown theoretically and numerically to be stable and convergent, with an optimal convergence rate.

A relaxed maximum-based normal $S$-iteration method (RMNSI) is proposed for solving generalized absolute value equations (GAVEs). The method combines a maximum-based fixed-point formulation with a constant relaxation parameter and avoids the selection of an auxiliary matrix. Assuming that $A+B$ is non-singular, both stages require linear systems with the same coefficient matrix $A+B$, allowing a single factorization to be reused throughout the iteration. A single spectral condition is established that guarantees unique solvability of the GAVE and global convergence from an arbitrary initial vector. An admissible range of the relaxation parameter is characterized, and $R$-linear convergence and error estimates are derived. Since the global condition can be conservative, a local convergence analysis based on the solution sign pattern is also presented. Numerical experiments on complementarity-derived, dense mixed-sign, and asymmetric ridge-regression problems are presented, and comparisons with several recent methods demonstrate the competitive efficiency and accuracy of RMNSI. The results further indicate that the preferred relaxation parameter depends mainly on the matrix structure and diagonal shift, while its dependence on the problem dimension is generally weak.

Alternating projections (AP) near a clean intersection of smooth manifolds can converge slowly when the Friedrichs angle is small. We develop the Fixed-Point Manifold Acceleration (FPMA) framework for {local convergence and polynomial acceleration of} smooth maps that fix a manifold pointwise. {Under suitable regularity and uniform normal spectral stability, FPMA establishes local R-linear convergence and a smooth limiting map inducing a local retraction. At fixed points, consistent polynomial acceleration with $p(1)=1$ preserves the selected limit's first-order response to normal initial perturbations whenever the transformed normal blocks are uniformly spectrally stable. Under sufficient smoothness, we give a curvature criterion for second-order induced retractions and show that stable consistent polynomial acceleration preserves both the quadratic retraction expansion and this criterion. {As a special case, we establish new local R-linear convergence guarantees for fixed-degree Chebyshev-accelerated exact AP and second-order induced retractions. For uniform Friedrichs angle $\theta_F\in(0,\pi/2)$, acceleration strictly improves the local convergence rate per AP evaluation. The large-degree limiting spectral factor at fixed $\theta_F$ predicts $O(\theta_F^{-1})$ evaluations for fixed relative error reduction as $\theta_F\downarrow0$, versus $O(\theta_F^{-2})$ for ordinary AP.}} {The framework also covers} inexact, relaxed, and generalized AP {under their respective regularity and stability conditions.}

The Independent Chip Model (ICM) converts the chip stacks of the players remaining in a poker tournament into finishing-place probabilities and prize equities. It is the standard model in tournament solvers, but its definition considers all possible finishing orders, and exact computation has been regarded as intractable for large fields. The same mathematics appears in other fields; for instance, the ICM is a Plackett-Luce ranking model with stacks as weights. We present DE-ICM, a deterministic algorithm that computes the placement probabilities of all $n$ players for all paid places in $O(M n^{2})$ time, where the number $M$ of quadrature nodes is fixed for the admissible inputs and the target accuracy. Its numerical scheme allows the target to be set close to the accuracy of double-precision arithmetic. The algorithm evaluates a classical integral representation in which, conditional on a player's exponential clock, the number of players ahead is Poisson-binomial. A double-exponential quadrature on a single grid evaluates every integral; its truncation errors have closed-form bounds, and the step size, which controls the discretization error, follows from the field size. A numerically stable deconvolution then recovers each player's leave-one-out coefficients in $O(n)$ time from one shared product. On exactly solvable instances with up to 4,000 players, the worst observed relative error of an equity is $4.0 \times 10^{-14}$ and the worst absolute error of a placement probability $1.1 \times 10^{-14}$. The full placement matrix for 1,000 players takes 0.2 seconds on one CPU core.

The goal of this article is to study optimal high-order Taylor-based polynomial quasi-Trefftz spaces for first-order time-harmonic Maxwell's equations. Compared with standard polynomial spaces, the purpose of quasi-Trefftz spaces is to retain the same high-order best approximation property but more efficiently: their dimension is considerably smaller. To achieve this goal, they approximate only smooth solutions to the governing equation rather than to approximate general smooth functions. While there is a natural definition of Taylor-based polynomial quasi-Trefftz spaces for scalar equations, the situation is fundamentally different for vector-valued equations, as the naive quasi-Trefftz property does not systematically define the smallest space with the desired best approximation property. This article identifies the optimal quasi-Trefftz spaces for first-order time-harmonic Maxwell's equations, proposes two approaches to construct quasi-Trefftz bases, and proves the spaces' best approximation property. The theoretical results, which strongly leverage polynomial Helmholtz decompositions and the exactness of the polynomial de Rham sequences, are validated with numerical experiments.

Focusing on scalar conservation laws, a TVD theory is developed for the fully discrete, single-stage Active Flux method. The classical choice of degrees of freedom is considered, yielding a third-order method. Starting from an adaptation of the concept of (discrete) total variation to Active Flux, a corresponding adaptation of Harten's lemma is achieved. Finally, the high-order method is endowed with parameters chosen such that they can be efficiently computed to fulfill the conditions of the lemma.

The minimum solid angle of a simplex is shown to be a shape measure in any dimension. While this has been known in dimensions two and three, the higher-dimensional generalization is new. We prove inequalities between the minimum solid angle of a non-degenerate simplex and its fullness. Up to dimension-dependent constants, the fullness is bounded from above by the minimum solid angle and from below by the square of the minimum solid angle. The exponents are optimal for three and higher dimensions. This work therefore generalizes the minimum angle condition. We also review several simplex shape measures and their relationships.

We study kernel functions induced by a linear time-invariant (LTI) stochastic differential equation (SDE), an approach that embeds prior knowledge about the function to be estimated directly into the state-space description of the SDE. This kernel-design perspective is relevant across the many settings where kernels play a central role, including kernel methods in statistics, Gaussian process (GP) regression in machine learning, and kernel-based regularization methods in system identification. The resulting class of SDE-induced kernels includes several kernels of practical interest, such as Matérn kernels with half-integer smoothness, spline kernels, and related kernels used in system identification. In contrast to existing linear-time GP regression methods for these kernels, which are typically based on Kalman filtering and Rauch–Tung–Striebel smoothing on the latent state, we work directly with the SSS representation of the observed covariance matrix. This yields linear-complexity recursive algorithms for the matrix operations required in GP regression, including factorization, linear solves, log determinants, posterior covariance, log-marginal-likelihood gradients, and posterior sampling; the same generator-level recursions extend, with only minor modifications, to the analogous computations in kernel-based regularization methods for system identification. Compared with previous uses of SSS representations in GP regression, we derive efficient recursive algorithms for analytic marginal-likelihood gradients in three application-relevant cases, and we develop new recursions for posterior computations and Matheron-type sampling. Numerical experiments verify the accuracy of the proposed recursions and illustrate posterior inference and sampling for representative SDE-induced kernels.