Papers

We study continuous-time and possibly high-dimensional stochastic control problems where drift coefficients and running reward functions are unknown. Due to these missing model primitives, we take the exploratory, reinforcement learning (RL) framework of Wang, Zariphopoulou, and Zhou(2020) with relaxed controls and entropy regularization. The objective is to develop theoretically grounded, efficient and scalable RL algorithms to learn both the optimal value functions (which also solve the exploratory HJB equation) and optimal exploratory feedback control policies. When the diffusion coefficients do not contain control, we employ probabilistic representations of both the optimal value function and its gradient based on an auxiliary state process depending only on the diffusion part of the original dynamics. With a delicate analysis on some properly defined mappings and their fixed points, this leads to the introduction of our policy iteration algorithms and their convergence. We demonstrate the performance of our algorithms through various numerical examples. Finally, we study a special control-dependent diffusion case where probability representation of the Hessian is called for.

Simulating elastic wave propagation in heterogeneous media presents significant mathematical and computational challenges, since resolving fine-scale material variations requires extremely fine spatial discretizations, while mixed stress–displacement formulations typically lead to large saddle-point systems that are not well suited for efficient time integration. Moreover, conventional multiscale reductions generally produce nontrivial coarse mass matrices, thereby requiring additional linear solves at every time step and diminishing the advantages of explicit schemes. To address these difficulties, we develop a time explicit multiscale method based on a multipoint stress control volume discretization. The central innovation is a unified spatial–temporal reduction strategy that simultaneously removes two dominant algebraic bottlenecks: the complex global fine-scale mixed solve and the repeated coarse-scale mass solve arising in conventional multiscale time stepping. More specifically, the stress and auxiliary rotation variables are eliminated locally, leading to a symmetric positive definite reduced stiffness operator, upon which a multiscale space is constructed using a density-weighted local spectral problem. The associated projected mass bilinear form preserves the physical mass inner product and gives an identity mass matrix under a density-orthonormal auxiliary basis, leading naturally to an explicit central-difference scheme. We establish discrete energy stability under a coarse-scale CFL condition and derive convergence estimates for both the density-weighted displacement and the locally recovered stress. Numerical experiments in heterogeneous media validate the theoretical results and demonstrate the effectiveness of the proposed method for elastic wave propagation.

We study boundary transmission and excursion geometry in the planar Lorentz mirror model through finite-volume numerical experiments. Cylinder penetration quantiles grow approximately linearly with circumference over the simulated range, with a systematic upward drift in the normalized quantiles when vacancies are present. Estimated slab-transmission probabilities yield exponential confinement bounds with quantified statistical confidence at the tested widths. Planar excursions conditioned to reach a distant boundary exhibit transverse ranges proportional to the exit radius, while their repeated-visit fractions vary strongly with mirror density. Endpoint-label collision probabilities substantially exceed a general gluing lower bound, although they account for only part of the observed same-column crossing probability. Finally, rerouting at repeated vertices preserves the exit while substantially shortening the sampled paths. These findings characterize finite-volume transport and excursion geometry without determining an asymptotic scaling exponent or resolving infinite-plane localization.

Any objects in space may encounter singularity due to the inverse-square law of the gravitational force, hindering the analysis and propagation of particle motion in space. Regularization in celestial mechanics removes such gravitational singularities in equations of motion by transforming the time and coordinates. This paper addresses a class of regularization methods known as global regularization for restricted three-body problems (R3BPs), which simultaneously removes the two singularities (e.g., Earth and Moon in cislunar space). It enables analysis, propagation, and prediction of particle motion throughout the domain of R3BPs, which is in stark contrast to regularizing each singularity individually. In particular, this paper develops a formulation for global regularization of the spatial elliptic R3BP (ER3BP) under external perturbations. The presented formulation leverages mathematical tools from quaternionic (hyper-complex) algebra and Hamiltonian mechanics. This approach enables us to (i) elegantly handle the four-dimensional coordinates arising from the dimension-lifting operation for spatial regularization, (ii) systematically derive the coordinate-momentum transformation while preserving the canonical structure and incorporating perturbing forces, and (iii) provide geometric insights for analyzing the regularity of the transformation, resulting Hamiltonian, and dynamical flow. Numerical examples demonstrate the validity of the mathematical results as well as significant numerical benefits of the regularization formulation in perturbed Earth-Moon ER3BP.

The generalized source integral equation (GSIE)approach is extended to enhance the compressibility of moment matrices for essentially-convex homogeneous penetrable scatterers. In the proposed indirect surface integral equation (SIE) formulation, the scattered field representations use a single distribution of generalized sources that radiate only weakly into the scatterer. A separate distribution of conventional sources is used for representing the interior fields. The wave interactions between surface subdomains via the generalized kernel are of a reduced effective dimensionality. The corresponding exterior integral operator moment matrix blocks exhibit enhanced rank deficiency. For lossy scatterers, this secures the low-rank (LR) compressibility of hierarchical matrix structures for the full system. For scatterers in a much slower background medium, significant savings with very delayed asymptotic scaling of the costs can be achieved. The GSIE is formulated for the transverse magnetic (TM) scattering problem. It is validated for representative examples and its superior compressibility is demonstrated.

We propose a localized core–tail Fourier–Laguerre frame method for approximation and model problems on the real line. The domain is decomposed into a finite core and two semi-infinite tails. Local Fourier extension is used in the core to resolve nonperiodic, oscillatory, and locally nonsmooth structures, while modulated Laguerre frames are used in the tails so that the Laguerre functions approximate only slowly varying decaying envelopes. The local Fourier extension component also provides two data-driven mechanisms: coefficient-energy indicators for internal edge detection and local frequency indicators for selecting tail modulation centers. We derive error estimates that separate the core approximation error, the modulated Laguerre envelope error, the effect of frequency mismatch, and the finite-tail truncation error. The analysis shows that the tail complexity is governed mainly by the residual phase after modulation rather than by the original carrier frequency. Numerical experiments demonstrate high accuracy for oscillatory, multi-frequency, and derivative-discontinuous functions, and a decaying model problem illustrates that the representation can be combined with differential operators through exact interface constraints.

Given a real square matrix $A$ and a nonempty closed target set of complex matrices $\mathcal E$, invariant by complex conjugation and containing real matrices, does the subset of $\mathcal{E}$ consisting of the matrices nearest to $A$ in the Frobenius distance always contain a real matrix? We give negative answers when $\mathcal{E}$ is either the set of normal matrices (solving a question posed by N. Higham) or the set of matrices whose eigenvalues lie in the closed left half-plane (solving a question posed by the first author and F. Poloni). For the latter problem, we also argue that the ratio between the real and complex distances is unbounded in every dimension $n\geq 3$, and this holds for every distance induced by a unitarily invariant norm. For both problems and $n \geq 3$, we show that a real input $A$ uniformly drawn from the unit sphere $\|A\|_F=1$ has no real minimizer over $\mathcal{E}$ with probability strictly between $0$ and $1$. The paper is complemented by some further results that are valid for more general target sets $\mathcal{E}$.

This paper introduces an online evolve–filter–relax (EFR) strategy for long-term stability of Operator Inference (OpInf) reduced-order models (ROMs) of complex fluid flow simulations. The main novelty of the new EFR-OpInf strategy is the use of online (i.e., at the learned ROM online evaluation level) spatial filtering inspired from large eddy simulation to significantly improve long-term stability and predictive performance of standard OpInf. Furthermore, the EFR-OpInf strategy reduces, and in some cases even eliminates, the need for standard $L^2$ regularization, while providing physical interpretability for the OpInf hyperparameters. The EFR-OpInf framework is modular, readily integrated into existing OpInf workflows, and accommodates a user-selected ROM filtering strategy. We demonstrate the EFR-OpInf's effectiveness using a fully non-intrusive projection-based ROM filter, a ROM differential filter, and a hybrid projection-differential ROM filter. The new EFR-OpInf models are evaluated on a high-Péclet-number convection–diffusion–reaction problem that embeds the variation in one parameter and two unsteady Navier–Stokes problems that focus on predictions beyond a training horizon: a transitional two-dimensional flow past a cylinder and a three-dimensional turbulent minimal channel flow. Across these three scenarios, EFR-OpInf can reduce prediction errors by up to an order of magnitude relative to standard OpInf. Moreover, EFR-OpInf remains stable over long prediction horizons in cases where standard OpInf diverges. Depending on the ROM filter used, EFR-OpInf's computational cost is comparable to that of standard OpInf.

We study torsional rigidity as a function of the labeled vertices of a convex polygon. Starting from the distributed second shape derivative, we derive the Hessian with respect to vertex coordinates. At a regular polygon, dihedral symmetry makes this matrix block circulant in radial-tangential coordinates, reducing its spectrum to the eigenvalues of Hermitian matrices of order two. We also derive an exact second-variation Galerkin identity and guaranteed functional residual majorants. Finite elements approximate the PDE solutions entering the Hessian, and FLINT/Arb provides the interval arithmetic needed for certification. In the scale-invariant setting, we certify exactly four zero eigenvalues generated by similarities and $2n-4$ strictly negative eigenvalues for $5\leq n\leq25$. The regular polygons in this range are therefore strict local maximizers, modulo similarities, of torsional rigidity divided by area squared. We prove broken regularity of the first and second material derivatives at each fixed regular polygon, and deduce an $O(h^2)$ error bound for the exact Galerkin Hessian and its eigenvalues on meshes fitted to the coarse fan.

A structure-preserving symmetric Galerkin boundary element method (SGBEM) is developed for the anisotropic multi-phase Mullins-Sekerka problem in $\mathbb R^2$. The bulk harmonic problems are eliminated through a discrete Dirichlet-to-Neumann operator assembled from region-wise SGBEM discretization, leading to a boundary-only formulation coupled to the evolution of a polygonal curve network with triple junctions. The Galerkin approximation of the resulting boundary problem is shown to be uniquely solvable, and the fully discrete scheme is unconditionally stable with respect to the time step as long as the evolving polygonal network remains in an admissible class. The accuracy of the proposed method is demonstrated by a convergence test under coupled space–time refinement against a similarity-rescaled exact solution consisting of three concentric circles.

We study two-point fractional boundary value problems with uncertain input data, where randomness may enter through the coefficients and boundary conditions. To quantify the resulting uncertainty in the solution, we employ the generalized polynomial chaos (gPC) framework and develop a stochastic Galerkin formulation of the problem. A particular focus of this work is the convergence analysis of the resulting approximation. Rather than imposing assumptions directly on the stochastic coefficients appearing in the gPC representation, we introduce minimal regularity assumptions on the input data and use them to establish the properties required for the convergence analysis. Based on these results, we prove the convergence of the stochastic Galerkin approximation to the corresponding gPC solution. Numerical experiments are presented to illustrate the theoretical findings and to investigate the influence of random coefficients and boundary conditions on the statistical behavior of the solution.

We derive localized maximum-norm bounds for bending moments computed by the Hellan-Herrmann-Johnson method for the clamped Kirchhoff plate problem. The main difficulty is to localize the discrete stress without leaving the HHJ space or violating its kernel constraint. Using symmetric-curl potentials, we construct a kernel-preserving localization and connect local Green stresses with a global discrete Green stress. This argument separates the local interpolation error from a weaker global pollution term. Under explicit regularity assumptions on the auxiliary problems, the resulting estimates give optimal pointwise convergence for bending-moment values and elementwise first derivatives. No logarithmic loss occurs for positive polynomial degrees, whereas the lowest-order value estimate retains a logarithmic factor. Under additional reflection symmetry, symmetric recovery improves bending-moment values for even polynomial degrees and first derivatives for odd degrees. The proved gains are one third and one half of an order, respectively. Numerical experiments confirm this parity dependence and exhibit gains close to one full order, exceeding those established theoretically.

This paper studies the strong convergence order of structure-preserving stochastic theta Milstein methods for a class of index-$1$ stochastic differential algebraic equations (SDAEs) with time-dependent singular matrices and non-globally Lipschitz coefficients. The singular matrix is allowed to vary in time while preserving a fixed differential algebraic splitting, and the drift and diffusion coefficients may exhibit superlinear growth. By exploiting the index-$1$ algebraic-differential decomposition of the exact solution, we identify the Milstein coefficient of the reduced stochastic differential equation directly in the original SDAE variables and establish the well-posedness and constraint preserving property of the proposed method for $\theta\in[1/2,1]$. Under a coupled monotonicity condition and suitable polynomial regularity assumptions, the method is proved to preserve the algebraic constraints at all time levels and to converge with strong order one in the root mean square norm. Numerical experiments confirm the structure-preserving property and the theoretical convergence order.

We develop and analyze a class of Hermite-type methods for linear advection-diffusion equations in the case of one-dimensional periodic domains. The methods evolve both nodal values and cell averages and are therefore referred to as hybrid-variable (HV) methods. Using Hermite interpolation theory, we construct the HV methods to arbitrary order of accuracy, and prove that the method is supraconvergent in the sense that the spatial order of accuracy can be up to two orders larger than the local truncation error. In the second part, we prove all central HV methods for linear diffusion equations and linear advection-diffusion equations are stable at the semi-discretized level. Next, we prove the convergence of all stable HV schemes and provide error estimate in the discrete $L^2$-norm of the semi-discretized solution, which extends the supraconvergence from a single mode to all sufficiently smooth initial data. The theoretical predictions on the accuracy and stability of the methods are verified by extensive numerical tests.

A tensor randomized average block coordinate descent method with heavy-ball momentum is proposed for the tensor least squares problem with respect to the t-product. Theoretical analysis for tensor block coordinate descent method is established and average techniques are applied to avoid the computation of tensor Moore–Penrose inverse. To further accelerate convergence, a heavy-ball momentum scheme is incorporated into the block coordinate descent method, where the step size and momentum parameters are adaptively determined by a two-dimensional minimal residual projection. Theoretical analyses give the convergence of the new method and provide an improved bound on the linear convergence rate. Numerical experiments further verify the efficiency of the proposed method in terms of the number of iterations and the better video recovery performance.

This paper introduces a unified three-phase data-driven twin framework for finite-horizon forecasting of physical quantities from coupled experimental measurements. The framework combines a new Hankel–Koopman finite-horizon energy decomposition with orthogonal modes and an inverse-calibrated multi-output nonlinear autoregressive model with exogenous inputs. The nonlinear dynamics are identified and simulated in the Hankel coefficient space over a prescribed calibration window. Candidate models are then evaluated after anti-diagonal recovery and reconstruction of the corresponding trajectories in the physical measurement space. The resulting simulated channel trajectories serve as exogenous inputs to a recursive forecasting model for the quantity of interest. In this way, the proposed framework integrates reduced-order representation, nonlinear dynamical identification, measurement reconstruction, and explainable forecasting. The mathematical analysis specializes Koopman delay-coordinate theory to the serialized multichannel Hankel lift and establishes its compatibility with the shifted Hankel representation. It also proves the orthogonality of the modes and the finite-horizon modal energy decomposition, together with the existence, uniqueness, inverse stability, and forward stability of the identified models under explicitly stated hypotheses. In a solar-power-plant case study, the framework yields consistently high correlations and low relative errors across all considered forecast horizons, demonstrating its ability to preserve the temporal evolution of the quantity of interest during recursive forecasting.

We develop a graph-compatible batched wavefront formulation for upwind discontinuous Galerkin discretisations of linear transport. For sweepable discretisations, the cell equations form a block lower-triangular system under a topological ordering of the directed upwind dependency graph. This causal structure can be shared across families of fixed-source problems whose local operators, sources and inflow data vary. We exploit this observation by grouping channels with a common dependency graph into graph-compatible sweep classes and carrying samples, right-hand sides, energy groups and sign-compatible angular ordinates as tensor dimensions within a common wavefront schedule. We prove that the resulting batched wavefront algorithm is algebraically equivalent to independent classical DG sweeps for every channel in a class, while exposing parallelism simultaneously across wavefront cells and channel dimensions. We realise the formulation as a GPU tensor program based on batched cell-local solves and quantify its work, depth and storage requirements, including the throughput–memory tradeoff introduced by sample microbatching. We further characterise the discrete parameter-to-observable map on fixed face-sign branches and show that reverse-mode differentiation through the wavefront algorithm reproduces the corresponding discrete adjoint action. Numerical experiments verify the expected DG convergence and adjoint consistency, demonstrate substantial GPU execution gains from sample batching, and illustrate the method in a material-shadowing uncertainty-quantification problem and a coupled multigroup C5G7/KAIST-inspired power-iteration workload.

Let $A$ be an $n\times n$ symmetric matrix with eigenvalues $\lambda_1\geq\cdots\geq\lambda_n$. Let $\tilde A:=A+E$, where $E$ is a symmetric noise matrix, and denote the eigenvalues of $\tilde A$ by $\tilde\lambda_1\geq\cdots\geq\tilde\lambda_n$. Bounding the perturbation $|\tilde \lambda_i -\lambda_i| $ is a central problem in linear algebra and numerical analysis. In this paper, we prove new perturbation results by exploring the actual interaction between $E$ and the eigenvectors of $A$. In the setting where $E$ does not act adversarially with respect to these vectors (for instance, if $E$ is random), we obtain a considerable improvement over Weyl's inequality. One can routinely extend these results to the Hermitian and rectangular settings. We obtain the new bounds as corollaries of a regional stability result, which provides a sufficient condition for a region on the real line to be stable (containing the same number of eigenvalues) after the perturbation. This result is of independent interest. We prove our regional stability result using contour integral analysis. Our main new technical ingredient here is the double-jump argument, which is robust and could be useful in many other situations involving Neumann series.

We derive computable upper bounds on the maximum wave speed in the Riemann problem for the compressible Euler equations with a composite equation of state, in which the pressure is the sum of a convex hydrodynamical part and a barotropic correction. The correction may destroy the convexity of the equation of state, so that the Riemann solution can contain composite shock-rarefaction waves. The bounds are obtained by comparison with the auxiliary Riemann problem for the hydrodynamical equation of state alone. The bounds hold irrespective of whether the outer waves are elementary or composite and require no case distinction on the wave type.

We develop controlled B-series and controlled Runge-Kutta methods as extensions of classical B-series and Runge-Kutta theory to controlled-driven rough differential equations. The driving signal is a path controlled by an underlying $\alpha$-Hölder step-$2$ rough path, with $\alpha\in(1/3,1/2]$, as arises when the output of one rough system drives another. Controlled B-series retain the classical rooted-tree and elementary-differential structure, while incorporating driver-dependent coefficients defined recursively through controlled rough integration. We establish finite-order expansions of the exact solution and its controlled Runge-Kutta approximation, and derive corresponding tree-based order conditions. Under suitable smoothness, solvability, stability, and boundedness assumptions, matching through tree order $p$ yields local error $O(|t-s|^{(p+1)\alpha})$ and global error $O(h^{(p+1)\alpha-1})$, provided $(p+1)\alpha>1$. The standard rough differential equation formulation is recovered when the controlled driver is the reference path itself, and the classical time-driven case recovers the usual B-series and Runge-Kutta methods. We further construct a simplified method using only increments of the controlled driver. Segmentwise canonical lifting of its piecewise linear interpolation gives explicit tree coefficients and reduces the third-order conditions to the classical ones. Assuming a Wong-Zakai solution-approximation rate $r_0>0$, the simplified method has global error $O(h^{\min\{r_0,(p+1)\alpha-1\}})$. Numerical experiments with fractional Brownian controlled drivers illustrate the resulting convergence behaviour.

Quantum computing offers the potential to solve computational mechanics problems with a more favourable complexity scaling than classically possible. To realise this potential, a quantum-compatible formulation of the finite element method is essential. We present such a formulation for Poisson problems with spatially varying coefficients on general domains discretised by an immersed uniform Cartesian grid. The resulting banded sparse system matrices are encoded as the sum of diagonal matrices with polynomial entries and their products with shift matrices. The sum of the matrices is formed using the linear combination of unitaries (LCU) technique. The resulting block-encoded system matrix can be implemented using a polylogarithmic number of elementary gates. For immersed boundaries, the block-encoded system matrix is further processed using a diagonal indicator matrix with zeros and ones. The indicator matrix is defined in terms of an integer-valued level-set function composed of simple geometric primitives via Boolean set operations. The required integer operations are computed on the fly using quantum arithmetic. The resulting block-encoded linear system of equations is solved using quantum singular value transformation (QSVT). The overall gate count, including the QSVT, is $O(N_{\text{tot}}^{2/d} \operatorname{polylog}(N_{\text{tot}}))$, where $d$ is the problem dimension and $N_{\text{tot}}$ the number of grid points. In numerical examples that confirm the theoretical complexity estimates, the solutions converge approximately linearly for immersed problems and quadratically for simply connected domains.

A matrix cascade is an iterative refinement process in which vector-valued data are repeatedly filtered by matrix-valued masks across successively finer scales. Such cascades arise in multichannel subdivision, Hermite refinement, multiresolution synthesis, and exponential-spline constructions, and in nonstationary schemes the masks may vary with scale. We develop a chronological stability framework for this setting, allowing both the transition operators and the associated matching spaces to change from level to level. If the chronological radius satisfies $\rho_{\rm ch}<1$ and the matching defects are summable, then the cascade converges uniformly and, for every $q\in(\rho_{\rm ch},1)$, \[ \|F_n-\Phi\|_\infty \lesssim q^n+\sum_{j=1}^n\delta_j q^{\,n-j} +\sum_{j>n}\delta_j . \] The time-weighted defect term records the level at which each perturbation is inserted, and the matching spaces need not converge; they may even alternate indefinitely. For compact piecewise-linear seeds, every finite ordered matrix cascade admits an exact fixed-width ReLU realization with $O(n)$ depth and parameters, and this linear depth order is worst-case optimal at fixed width. Consequently, exponentially decaying defects yield certified approximants of size $O(\log\varepsilon^{-1})$. We further construct stable two-channel classes in which the chronological radius and the leading metric-entropy coefficient are independent invariants, showing that stability does not determine finite-accuracy description complexity. The framework therefore provides tolerance-controlled finite representations for nonstationary multiresolution and exponential-spline refinement.

This work studies subspace embeddings obtained by two normalized real Walsh transforms, two independent sign diagonals, and uniform coordinate sampling without replacement. The main result shows that the prescribed sample size $k=\min\{n,\lceil Cr/\varepsilon^2\rceil\}$, for a universal constant $C$, suffices to preserve all squared norms on each fixed $r$-dimensional subspace within $1\pm\varepsilon$ with probability at least $0.99$. The result holds for every ambient Walsh dimension and all ranks, and answers Problem TR-01 in the Open Problems in Numerical Linear Algebra repository. The proof controls joint entry cumulants of the transformed projection through connected graph contractions and Walsh character identities. These estimates then bound the expected trace of even powers of a product of centered projections. A two-projection decomposition converts this estimate into control of both spectral edges. Bernoulli sampling at arbitrary densities and a deterministic upper bound near full sampling yield the prescribed number of coordinates.

This paper introduces adaptive moving mesh methods for the numerical simulation of axisymmetric mean curvature flow, addressing both isotropic and anisotropic cases. The methods are developed within the framework of the mesh equidistribution principle, where a carefully designed tangential velocity is employed to dynamically redistribute mesh points during the evolution. To accurately capture the key geometric features of the evolving interfaces, we select monitor functions based on the curvature $\kappa$, its arc-length derivative $\kappa_s$, and the squared curvature $\kappa^2$. These monitor functions can be flexibly tailored to suit different problem settings and play a vital role in determining the resulting mesh quality and numerical accuracy. Spatial discretization is performed using central finite differences, while temporal integration is handled with first- and second-order time-stepping schemes, including the BDFk ($k=1,2$) and Crank-Nicolson methods. Additionally, a Lagrange multiplier approach is incorporated into the adaptive system to enforce the underlying geometric constraint, resulting in energy-stable numerical schemes. Numerical experiments confirm the convergence and energy stability of the proposed methods. More importantly, the results clearly show that the proposed methods offer significant advantages in complex geometric evolutions: by utilizing appropriately designed monitor functions, the adaptive methods achieve dynamic redistribution of mesh points, efficiently capturing localized geometric features, significantly improving numerical accuracy, and effectively preventing mesh degeneration, particularly in anisotropic cases.

We propose and analyze a non-iterative domain decomposition time integrator for the linear acoustic wave equation. The method is based on an overlapping decomposition in space and uses a localized implicit prediction to provide artificial boundary data at the interfaces. We use conforming finite elements of arbitrary order to discretize in space and allow spatially varying material coefficients. We prove an error bound relative to the global Crank–Nicolson approximation that decays exponentially with the size of both the subdomain and prediction overlaps. Compared to earlier work, this new method avoids a CFL-type step size restriction and allows the use of finite elements of arbitrary order. The analysis combines a discrete Saint-Venant principle with new decay and truncation estimates for the localized prediction. Numerical experiments support the analysis and demonstrate stable behavior for all tested time step sizes, accuracy comparable to that of the global Crank–Nicolson approximation, and the potential for efficient parallel execution.

Operator learning is formulated on function spaces, but training data are typically available only through finite-dimensional representations. In encoder–decoder architectures, a matrix-valued kernel on the encoded space induces an operator-valued kernel on the original function spaces, and the corresponding reproducing kernel Hilbert spaces are isometrically isomorphic. As the input and output resolutions increase, the induced kernels converge to a limiting kernel, in the sense of operator-norm convergence of their associated integral operators, allowing regularity assumptions to be stated independently of the encoding resolution. For regularized stochastic gradient descent, we establish upper bounds for decreasing and fixed step sizes, separating the encoding and regularization terms from optimization terms of order \(t^{-\theta}\) and \(T^{-\theta'}\), respectively, for any \(\theta,\theta'\in(0,1)\). We further prove lower bounds showing that these encoding-induced terms are generally unavoidable. The analysis is further extended to encoder–decoder neural networks through the limiting neural tangent kernel (NTK), yielding error bounds with an additional finite-width term and polynomial parameter and sample complexity guarantees when the encoding errors decay algebraically. The framework covers matrix-valued kernels constructed from radial and dot product kernels, NTKs arising from wide encoder–decoder neural networks, and encoder–decoder pairs based on Fourier, Legendre polynomial, wavelet, PCA, or pointwise sampling representations.

We study ground-state solutions of the 4th order nonlinear Schrödinger equation with mixed dispersion via the elliptic profile equation $\Delta^2 Q + 2a \Delta Q + bQ - |Q|^\alpha Q=0$, focusing on symmetry breaking of least-action ground states in 2D in the regime $a=1$, $b=1+\epsilon$, $0<\epsilon \ll 1$, where nonradiality was proved by Lenzmann and Weth [26]. We establish a quotient-to-mass relation converting small-$\epsilon$ asymptotics of the Weinstein quotient into mass asymptotics, in any dimension and symmetry class. In 1D we prove the sharp small-$\epsilon$ rate for the quotient for every $\alpha>0$ and deduce the rates for the mass, action and potential norm. In two and higher dimensions the same relation, combined with the quotient asymptotics of Lenzmann-Weth [26] and Mandel-Oliveira e Silva [27], gives the mass rates for the Knapp-type nonradial and the radial ground states. The resulting threshold $\alpha_K(d)=8/(d+1)$ coincides for $d\leq3$ with the existence of normalized minimizers threshold of Fernández-Jeanjean-Mandel-Mariş [16], connecting the least-action and mass-constrained formulations. We complement this with a detailed numerical study. For the cubic nonlinearity we construct a nonradial branch from a Knapp-type example, compare it with radial and angular-mode branches, and verify the predicted rates. For the quartic nonlinearity, turning points and branching are present in the mass-energy diagrams in radial or two- and four-peak nonradial solutions. Here, although least-action ground states are nonradial for small $\epsilon$, the normalized ground states of the same mass appear radial, a consequence of branching. For $\alpha\geq4$ we find only radial ground states, consistent with the conjectured sharp form of the Stein-Tomas inequality on $\mathbb S^1$, equivalent to the radial and nonradial masses having the same leading coefficient.

LP-Newton solves a linear program through a sequence of nearest-point problems. Given an interior feasible point for an inequality-form LP, we construct the compact hull formed by the origin and its normalized constraint rows. Polar LP-Newton (P-LPN) follows the objective ray to the boundary of this row-polar hull, where it recovers a primal-dual optimum or a recession direction proving unboundedness. For rational row-polar data, we prove an outer-iteration bound quadratic in dimension and linear in binary input length, excluding inner Wolfe work. Restarting Wolfe at every outer iteration discards its terminal corral even though the hull is unchanged and the next target lies on the same ray. Corral-repair P-LPN (CR-P-LPN) instead repairs that corral and uses it to start the next projection. Verification over the full hull preserves P-LPN's projections, outer targets, and LP conclusion in exact arithmetic. Experiments in Julia and MATLAB show where repair helps. After common initialization, median P-LPN times across 180 single-LP tests are 1.62 times the corresponding CR-P-LPN times in Julia and 1.55 times in MATLAB. On three problems with 5,000 to 100,000 rows, the corresponding geometric-mean factors are 2.03 and 1.46. A separate end-to-end study includes initialization on 19 application-derived LPs. CR-P-LPN has lower median total time than the faster cold HiGHS mode, either dual simplex or the interior-point method with crossover, on 18 of the 19 workloads in each language. The geometric means of the ratios of HiGHS time to CR-P-LPN time are 4.93 and 3.06. Component comparisons identify terminal-corral reuse as the useful component, whereas candidate-first ordering alone shows no consistent gain. Repair is not always beneficial: it is slower on Klee-Minty tests, and an author-implemented simplex is faster on the retained Netlib models.

For computing the ground state of the defocusing Gross–Pitaevskii energy, we propose and analyze two efficient schemes based on the Davis–Yin three-operator splitting, which treats the potential and interaction terms explicitly, and the kinetic energy by a resolvent. One iteration costs one inversion of $I-\gamma\Delta$ for the first scheme, and of $I-\gamma\Delta+\gamma V_1$ for the second, with $V_1$ denoting the separable part of the potential, and on structured meshes both operators can be inverted by simple fast GPU solvers. For monotone discrete Laplacians, including the second-order finite difference scheme and the lumped linear finite element method on simplicial meshes with suitable angle conditions, we prove global convergence to the unique positive discrete ground state, for every positive normalized initial vector, for any constant step size below an explicit threshold. In contrast, the methods previously proven to converge globally to the ground state all invert a more difficult elliptic operator that depends on the iteration variable. In three-dimensional tests with up to $999^3$ unknowns on one GPU, a simple variable step size rule makes the proposed splitting schemes efficient in practice, comparable in wall-clock time to Riemannian conjugate gradient methods that also invert only a shifted Laplacian operator, and much more robust with respect to the choice of the initial guess.

High-order entropy stable (ES) methods have attracted much interest in the computational fluid dynamics (CFD) research community over the past decade, as we seek to simulate practical engineering applications with high accuracy and robustness. As industry delves further into supersonic flight and space travel, we need schemes capable of modelling flows with multiple chemical species that interact and mix within the fluid. As such, there is a considerable amount of effort currently being focused on developing high-order ES schemes for multi-species (MS) flow. A key ingredient for these entropy stable schemes is the two-point flux, which determines the properties of the resulting high-order discretization. As such, the numerical scheme can possess important properties such as entropy conservation and kinetic energy preservation only if the two-point flux also mirrors these properties. At present, in the literature, the two-point fluxes used for high-order ES multi-species solvers are not kinetic energy preserving (KEP). In this short communication we develop a two-point flux that is both entropy conserving (EC) and kinetic energy preserving by adapting Ranocha's single-species methodology for multi-species. The newly proposed EC and KEP flux is tested against existing two-point fluxes in the literature using the density pulse advection and Taylor-Green vortex cases. It is shown that EC and KEP flux retains both properties while the other fluxes only retain one of the desired properties. Since the two-point flux dictates the properties of the resulting scheme, the proposed flux is the ideal option for simulating multi-species compressible flow.