Approximation

Physics-informed neural networks and finite element methods provide two different paradigms for the numerical approximation of partial differential equations: the former are commonly trained by minimizing pointwise strong residuals, whereas the latter are naturally built from weak variational formulations and the finite-dimensional systems obtained after discretization. In this work, we introduce a common framework based on the discretization of functional Gauss--Newton problems by finite families of linear measurements. We show that, through an appropriate duality pairing, the linear measurements can be represented by test functions. The resulting Gauss--Newton system is then precisely a Petrov--Galerkin discretization of the linearized functional problem. This perspective recovers pointwise collocation and natural-gradient constructions as particular cases, while making the choice of test functions an explicit algorithmic design choice. We specialize this framework to elliptic problems, where it naturally leads to weak residual formulations and to a hybrid finite element--neural construction acting on complementary approximation spaces. Numerical experiments support the proposed framework and demonstrate the effectiveness of weak Gauss--Newton formulations and hybrid finite element--neural approximations.

Recent studies have shown that smooth functions can be well approximated by ReLU neural networks with path norm constraint on the weights. We extend these results from uniform approximation to approximation in Sobolev norm. Specifically, we analyze how well Sobolev functions in $W^{n,p}$ can be approximated by neural networks with width $W$, depth $L$ and path norm bounded by $K$, when the approximation error is measured in the $W^{1,p}$-norm. For shallow networks with depth $L=1$, we derive the approximation error bound $\mathcal{O}(\max\{W^{-(n-1)/d}, K^{-(n-1)/(s-n)}\})$, when the smoothness index satisfies $n<s=(d+3)/2$ and the input is $d$-dimensional. For deep networks, we remove the restriction on the smoothness by showing that the approximation bound $\mathcal{O}(K^{-(n-1)/(d+d/p+1)})$ holds if the width $W$ and depth $L$ are sufficiently large.

Spectral co-clustering is a useful tool for discovering latent structure in word-document matrices, but its reliance on singular value decomposition (SVD) can make standard formulations expensive on high-dimensional data. This paper presents two randomized approximations for normalized spectral co-clustering of bipartite text data when the numbers of document and word clusters may differ. The first method uses randomized SVD through random projection, while the second combines partial SVD with element-wise random sampling. Across real-world and synthetic datasets, both methods reduce runtime relative to the full-SVD baseline, but their behavior depends on matrix sparsity. The random projection method is the more reliable approximation across the tested settings, whereas the sampling-based method is most useful on denser matrices and provides limited benefit on already sparse text data. These results show that randomized approximations for spectral co-clustering should be selected according to the underlying structure of the data.

Constructing Lyapunov functions for nonlinear dynamical systems is a central problem in stability analysis, yet remains challenging. Lyapunov functions are commonly characterized as solutions to first-order partial differential equations (PDEs), but these solutions are typically obtained for single systems, limiting their reuse across systems. In this paper, we study the Lyapunov solution operator that maps a vector field to the corresponding Lyapunov function defined by a dissipation-based Lyapunov PDE. We establish that, on compact subsets of the domain of attraction and under exponential stability assumptions, this operator is well-defined, unique, and continuous with respect to perturbations of both the vector field and the dissipation function. These results provide a theoretical foundation for approximating Lyapunov functions uniformly over families of nonlinear systems. Building on these theoretical foundations, we employ Fourier Neural Operators (FNOs) as a data-driven approximation of the Lyapunov solution operator. Numerical experiments demonstrate that a single trained operator can accurately approximate the numerical Lyapunov functions across parameterized families of dynamics. This illustrates the potential of neural operators for approximating Lyapunov functions.

Approximating the input-output behavior of a multivariable black-box function from limited data is challenging when blind to the importance of its inputs and their interactions. We introduce total sensitivity kernels (TSKs), a method based on families of weighted ANOVA kernels that learn and adapt to this multivariable structure. TSKs parameterize the weights on each multivariable component of the target function by factors for each input. We propose learning these factors directly from function evaluations by selecting the reproducing kernel Hilbert space (RKHS) in which the target function has minimum norm. Under suitable conditions, we show that this norm-minimization problem admits a unique solution, and we establish consistency of a finite-data formulation based on minimum-norm interpolation. The learned TSK factors characterize the participation of individual inputs across interactions and main effects, providing a kernel-dependent notion of input sensitivity related to total Sobol indices. Numerical experiments demonstrate that adapting the kernel to learned multivariable structure can substantially improve approximation accuracy over a standard product kernel.

We study the uniform approximation of smooth scalar-valued functionals on an infinite-dimensional separable Hilbert space by ReLU neural networks. A key feature in deep learning for functional data is the varying importance of different coordinates/dimensions. Representing the functional input in a basis expansion, we quantify the importance of each coordinate through both the magnitude of its corresponding basis score and the directional sensitivity of the target functional. Our analysis combines coordinate truncation, anisotropic partitioning, local Taylor approximation, and ReLU network realization, while allowing unrestricted interactions among the retained coordinates. We establish a general nonasymptotic upper bound for the uniform approximation error and a complementary pseudo-dimension-based lower bound for the worst-case approximation error. Under generalized exponential coordinate decay $w_ds_d\asymp\exp(-cd^ρ)$, with $ρ>0$, the upper and lower bounds match at the leading order, which is stretched-exponential in the logarithm of the network size budget, and thus yield the nearly optimal approximation rate. This is the first work to characterize neural network approximation error for infinite-dimensional functional inputs explicitly through the joint dimensional decay of coordinate magnitudes and directional sensitivities.

Multi-task learning (MTL) requires navigating unavoidable trade-offs among competing objectives. This paradigm is frequently formulated as multi-objective optimization (MOO), where the scalarization is favored to reduce an MOO problem to a single objective. We empirically find that existing merit-function-based scalarization approaches are sensitive to the relative scales of different objectives in practical MTL, where task losses commonly differ by orders of magnitude. The optimization process often favors objectives with larger scales even though the underlying Pareto optimal solutions remains invariant to rescaling (i.e., multiplying an objective by a positive constant). To address this issue, we propose Scale-Invariant Merit-function-based Scalarization (SIMS) for MTL. Specifically, SIMS adopts a transformation-induced merit function to convert the MOO problem of MTL to a single objective that renders optimization invariant to the magnitudes of losses. Theoretically, we prove that the requirement for scale invariance uniquely determines this transformation to be logarithmic. We further show that this general transformation-induced merit function preserves weak Pareto optimality and admits a smooth surrogate with controllable approximation error. Extensive experiments on representative multi-task benchmarks demonstrate that SIMS consistently outperforms existing scalarization methods and achieves state-of-the-art performance.

Stochastic gradient descent (SGD) admits diffusion approximations that replace the complicated randomness of stochastic gradients by Gaussian noise, providing a powerful tool for understanding its dynamics and long-time behavior. We investigate whether an analogous approximation principle holds for optimization over probability measures, where the objective is a functional defined on the Wasserstein space P2. The nonlinear geometry and infinite-dimensional nature of P2 prevent a direct extension of the classical Euclidean theory. Using Lions differentiability, we lift the problem to a linear Hilbert space, where higher-order differential calculus becomes available. We then construct a Gaussian random-field approximation whose velocity field matches the mean and covariance of the original stochastic gradient. By exploiting this moment matching through higher-order Taylor expansions, we show that the Gaussian approximation captures the SGD dynamics with second-order weak accuracy. Our result provides a rigorous foundation for replacing sample-driven randomness by analytically tractable Gaussian fluctuations in stochastic optimization over probability measures.

Multidimensional item response theory (MIRT) is a statistical test theory that precisely estimates multiple latent skills of learners from the responses in a test. Both compensatory and non-compensatory models have been proposed for MIRT: the former assumes that each skill can complement other skills, whereas the latter assumes they cannot. This non-compensatory assumption is convincing in many tests that measure multiple skills; therefore, applying non-compensatory models to such data is crucial for achieving unbiased and accurate estimation. In contrast to tests, latent skills will change over time in daily learning. To monitor the growth of skills, dynamical extensions of MIRT models have been investigated. However, most of them assumed compensatory models, and a model that can reproduce continuous latent states of skills under the non-compensatory assumption has not been proposed thus far. To enable accurate skill tracing under the non-compensatory assumption, we propose a dynamical extension of non-compensatory MIRT models by combining a linear dynamical system and a non-compensatory model. This results in a complicated posterior of skills, which we approximate with a Gaussian distribution by minimizing the Kullback-Leibler divergence between the approximated posterior and the true posterior. The learning algorithm for the model parameters is derived through Monte Carlo expectation maximization. Simulation studies verify that the proposed method is able to reproduce latent skills accurately, whereas the dynamical compensatory model suffers from significant underestimation errors. Furthermore, experiments on an actual data set demonstrate that our dynamical non-compensatory model can infer practical skill tracing and clarify differences in skill tracing between non-compensatory and compensatory models.

Consistent submodular maximization studies the tradeoff between solution quality and stability when elements arrive over time. For a monotone submodular objective, which models diminishing returns, an algorithm maintains a set of at most $k$ available elements and changes only $O(1)$ elements after each insertion. Dütting et al. [2025] established a tight $2/3$ approximation with unrestricted computation and a polynomial-time $0.51$ approximation. They left open at STOC 2025 whether efficient algorithms can match the offline $1-1/e$ guarantee. We resolve this problem by proving that the supremum approximation achievable with polynomially many value queries and worst-case constant recourse is $$β=2-\sqrt2\approx0.5858<1-1/e.$$ For every $\varepsilon>0$, our randomized algorithm attains $β-\varepsilon$ with $O(\varepsilon^{-2})$ changes per insertion. Any fixed improvement requires exponentially many queries before one critical insertion or linear recourse of $Ω(k)$ changes at that insertion, even with unlimited queries afterwards. This gap quantifies the cost of consistency: the current oracle hides which elements will be needed after an arrival. We also determine the exact curvature-dependent threshold $1-(\sqrt2-1)\vartheta$, attain $1-1/e-\varepsilon$ for weighted coverage with $O(\varepsilon^{-1})$ recourse, and separate the existence of universal future-price certificates from their efficient computation. Our algorithm has a bounded-bit polynomial-time implementation for polynomial-bit rational oracle answers; the lower bound uses only logarithmic-bit rational answers.

We characterize the sharp structure-agnostic minimax risk for coefficient estimation in the partial linear model when the outcome and treatment nuisances are learned by two distinct black-box learners, which resolves the open problem in double machine learning posed by Gu (2025). For each nuisance $q\in\{μ,π\}$, we characterize the available learner by an approximation-error budget $a_q$ and a stochastic-error budget $s_q$, with the latter controlled through localized Rademacher complexity. Writing $\mathcal E_n$ for the minimax mean-squared error, we show that $$\mathcal E_n\asymp1\wedge\left\{\frac1n+\left(a_μa_π+\min\left\{a_πs_μ+s_π^2,\,a_μs_π+s_μ^2\right\}\right)^2\right\}.$$ The main new ingredient is a novel lower bound for the general two-learner problem. Our proof constructs four finite-mixture testing experiments using orthogonal code functions. Across these experiments, the hidden perturbations are placed outside both learner classes, outside only the treatment learner class, outside only the outcome learner class, or inside both learner classes. These four configurations capture, respectively, the interaction between the two approximation errors, the two asymmetric interactions between one learner's approximation error and the other learner's learning error, and the joint estimation difficulty of learning both nuisances. Combining the four resulting lower bounds yields the displayed rate, which matches the latest upper bound in Gu (2026). Our result shows that standard double machine learning can overstate the intrinsic difficulty of target estimation and provides a target-specific principle for learner selection: approximation error and stochastic complexity must be jointly balanced across the two nuisance learners rather than optimized separately.

Fairness in clustering has attracted sustained research interest, motivated by the need to ensure equitable representation of protected groups in machine learning applications. We study fair $k$-means clustering in Euclidean space, where the proportion of each protected group in every cluster must lie within specified lower and upper bounds. These constraints make it challenging to determine both cluster centers and point assignments. We propose an approximation algorithm that combines a linear programming relaxation with geometric transformations of the input to construct candidate center sets. Given a $ρ$-approximate algorithm for weighted $k$-means and any $ε>0$, our algorithm returns a fractional solution whose cost is at most $1+(3-1/Γ)ρ+O(ε)$ times the optimal integral fair cost, where $Γ\approx6.357$ is an upper bound on the integrality gap of the standard Euclidean $k$-means LP. With a PTAS as the subroutine, the approximation ratio becomes $3.8427+O(ε)$, improving the previous factor of $5+O(ε)$ to below $4$. The solution satisfies all fairness constraints exactly and can be rounded to an integral assignment with a bounded additive violation of fairness and no increase in cost. The same approximation guarantee extends to the $k$-sparse Wasserstein barycenter problem.

Model compression is key to mitigate deployment challenges of ever growing machine learning models. In this area of research, singular value decomposition (SVD)-based compression offers a compelling trade-off between computational efficiency and model accuracy. Fisher-weighted SVD in particular provides principled, loss-aware compression. However, we find that improving the fidelity of Fisher approximation used in the compression is poorly predictive of post-compression accuracy for Vision Transformers (ViTs). Motivated by this observation, we propose FACTS, a structured Fisher Approximation tailored to Compressing ViTs with Fisher-weighted SVD, which enforces token-local aggregation while preserving within-token activation-gradient dependence. Additionally, we introduce a fast Constrained Rank Search (CoRS), that optimizes layer-wise rank allocation while adhering to a fixed floating point operation (FLOP) constraint. Extensive experiments across ViTs and hybrid architectures demonstrate that FACTS consistently improves accuracy-efficiency trade-offs without requiring finetuning. Notably, it outperforms the strongest SVD baseline by up to +5.8 percentage points (p.p.) Top-1 on Swin-B, with further gains driven by our search method. Code is available at https://github.com/MoritzTho/FACTS.

Rigorous results show that feedforward neural networks can overcome the curse of dimensionality in the numerical approximation of high-dimensional partial differential equations (PDEs), but comparatively little is known about residual neural networks (ResNets) in the nonlinear PDE setting. We prove that ResNets overcome the curse of dimensionality in the numerical approximation of solutions of semilinear heat equations with globally Lipschitz continuous, gradient-independent nonlinearities: under polynomial growth and network approximability hypotheses on the PDE data, there exist $η\in(0,\infty)$ and ResNets $Ψ_{d,\varepsilon}$, $d\in\mathbb{N}$, $\varepsilon\in(0,1]$, with at most $ηd^η\varepsilon^{-η}$ parameters whose realizations approximate the solution in dimension $d$ with an $L^2$-error of at most $\varepsilon$. The proof represents one deterministic realization of a multilevel Picard estimator by a ResNet whose shortcut connections transmit the spatial variable and a scalar accumulator, while the residual branches successively add the summands of the estimator. For ridge-sum initial conditions, admissible sigmoidal activations, and globally Lipschitz truncations of the nonlinearity, we obtain, for every $ξ>0$, the explicit bound $C_ξd^{4+ξ}\varepsilon^{-(3+ξ)}$ on the number of parameters.

Mixture-of-experts (MoE) architectures increase model capacity by combining a collection of expert predictors through input-dependent routing, while often activating only a small subset of experts for each input. Despite their growing importance in modern large-scale models, the statistical roles of their design choices, especially routing, sparse activation, and shared experts, remain only partially understood, as existing theory has largely focused on parametric or correctly specified MoE models. In this paper, we view MoE as a form of localized aggregation and show how this localization reshapes the approximation-estimation-computation tradeoff. We derive oracle risk bounds for learning dense and sparse routing with evolving experts, separating approximation, expert-learning, and router-estimation errors, and characterize how sparse Top-K routing can retain the benefits of localized aggregation while controlling per-input computation. We also interpret gating through the geometry of input space, relating routing performance to regions of local expert advantage, and show how shared experts, as adopted in architectures such as DeepSeekMoE, can extract common predictive structure so that routed experts focus on residual local variation. Together, these results provide a unified statistical framework for understanding MoE through input-dependent expert aggregation, in which expert specialization and computational tradeoffs are governed by local predictive structure.

This paper introduces the Sierpiński-Knopp (SK) Wasserstein distance, a fast metric between persistence diagrams. The SK-Wasserstein distance, denoted $d_{\mathrm{SK}}$, maps diagram points and their diagonal projections to the unit interval via the Sierpiński-Knopp space-filling curve on the upper diagonal triangle. The encoded point sets are then efficiently matched via one-dimensional optimal assignment, in $O(N\log N)$ steps, yielding an explicit diagonal-aware point assignment between the two input persistence diagrams. We show that the SK-Wasserstein distance controls the classical $2$-Wasserstein distance between diagrams, admits an explicit isometric embedding into a Hilbert space, and induces a positive-definite Gaussian kernel, making the resulting geometry directly compatible with Euclidean and kernel-based learning methods. A tighter surrogate dissimilarity, noted $W_Γ$, is also introduced based on the point assignments along the curve. Experiments on 12 scientific collections comprising 227 diagrams show median per-collection speedup of $d_{\mathrm{SK}}$ over state-of-the-art approximations of $W_2$ is $626\times$, while the aggregate speedup over the full benchmark is $2100\times$. Average-linkage partitions obtained from $d_{\mathrm{SK}}$ and $W_Γ$ each exactly match the corresponding $W_2$ partition on 8 of the 12 collections. Hilbert $k$-means and Gaussian spectral clustering, both based on $d_{\mathrm{SK}}$, achieve mean adjusted Rand indices (ARI) of $0.756$ and $0.800$, respectively, with respect to the benchmark reference partitions, compared to $0.750$ obtained by average linkage on $W_2$. The Gaussian $d_{\mathrm{SK}}$ kernel supports other kernel-based analysis tasks, as illustrated by its use for contiguous segmentation of ordered diagram collections in our experiments.

Phase retrieval in broadband coherent anti-Stokes Raman spectroscopy (BCARS) is an ill-posed inverse problem. The Raman-like signal is encoded in the imaginary part of the resonant susceptibility, which mixes coherently with a non-resonant background (NRB) that varies across acquisitions. We introduce an inverse physics-informed neural network (iPINN) that predicts Lorentzian peak parameters from raw BCARS spectra and reconstructs the resonant susceptibility through a differentiable analytical forward model. A transformer encoder assigns spectral features to 24 learnable peak slots, and a multi-view consistency loss enforces invariance across NRB pattern, NRB strength, and noise. Unlike direct spectral regression approaches, the method retains accuracy under varying acquisition conditions. On a public benchmark, iPINN achieves the lowest error among the tested baselines (MAE 0.016 vs. next-best 0.046). On 28 zero-shot test spectra acquired across seven solvents and four focal positions, accuracy is depth-invariant in five of seven solvents. These results show that inverse parametric prediction with a differentiable physical decoder supports robust phase retrieval across measurement conditions.

We prove that the Hypergraph Neural Network, an invariant architecture with 3-body message passing, is a universal approximator for potential energy surfaces. Our main contribution is a multi-layer completeness theory. We show that $L$ layers of message passing on sparse, cutoff-based graphs achieve the same representational power as having access to the full $L$-hop neighborhood, provided the configurations are generic, satisfy an overlap condition and a connectivity condition. This provides the first rigorous justification for the common practice of using multi-layer message passing with a per-layer cutoff smaller than the physical interaction range, the setting used by virtually all practical graph neural network based machine-learned interatomic potentials. As immediate consequences, we show that both DPA3 and CHGNet architectures inherit universal approximation.

Many parameter-efficient methods generate the parameters of a large neural network from a low-dimensional latent representation. Given an architecture $Φ$ with $P_Φ$ parameter slots, we write $\boldsymbolθ_f=\mathcal{G}(\boldsymbolξ_f)$, where $\mathcal{G}\colon\mathbb{R}^M\to\mathbb{R}^{P_Φ}$ is a parameter generator and $\boldsymbolξ_f\in\mathbb{R}^M$ is a latent representation of the target function $f$. The architecture $Φ$ and the generator $\mathcal{G}$ are shared across the entire target class, while each target $f$ is represented by its own latent vector $\boldsymbolξ_f$, with $Φ_{\mathcal{G}(\boldsymbolξ_f)}$ approximating $f$. This framework encompasses hypernetworks, low-dimensional parameterizations, parameter-efficient adaptation, and model compression. Understanding the tradeoff between the latent dimension $M$ and the network budget $P$ is therefore fundamental to characterizing the expressive efficiency of these methods. We study this tradeoff for affine generators and fully connected ReLU architectures. More precisely, optimizing jointly over architectures $Φ$ satisfying $P_Φ\leq P$ and affine generators $\mathcal{G}:\mathbb{R}^M\to \mathbb{R}^{P_Φ}$, we prove that the optimal worst-case uniform approximation error over the unit ball of $α$-Hölder functions on $[0,1]^d$, where $0<α\leq1$, has the sharp order $\bigl(P\min\{M,P\}\bigr)^{-α/d}.$ In particular, our result shows that even a fixed-dimensional latent space suffices to achieve vanishing approximation error as the network budget increases.

Bilevel mixed-integer linear optimization problems model hierarchical decision processes in which a leader anticipates the optimal response of a follower. Although expressive, these problems are computationally challenging because lower-level optimality is embedded in the leader's feasible region. Value-function reformulations replace the nested follower optimization with a constraint involving the follower's optimal value, but evaluating this value function exactly can itself be expensive. This paper introduces Graph4BiLO, a graph neural network (GNN) approach for learning bilevel value functions from variable--constraint graph representations. In contrast to fixed-length multilayer perceptron (MLP) representations, the GNN uses shared message-passing parameters and can therefore be applied across multiple problem sizes with a single trained model. The learned ReLU network is encoded exactly as mixed-integer linear constraints and embedded in an approximate single-level formulation. A repair step subsequently re-solves the follower problem for the selected leader decision to recover a bilevel-feasible follower response. We evaluate Graph4BiLO on knapsack interdiction instances with 20--100 items against the exact MibS solver and the learning-based Neur2BiLO method. Graph4BiLO obtains objective values comparable to Neur2BiLO across all tested sizes while avoiding size-specific neural networks. An additional out-of-distribution experiment demonstrates zero-shot transfer from 20-item training instances to previously unseen 40- and 60-item instances. However, embedding message passing at every graph node substantially increases the resulting mixed-integer formulation size and solve time. These results identify a central tradeoff between size-generalizable graph representations and the computational cost of embedding GNNs within optimization models.

Token prediction is a central pre-training objective for modern language models. Despite its empirical success, why token prediction learns broadly useful representations remains incompletely understood. We develop a statistical framework connecting token prediction with representation geometry, encoder approximation, and downstream performance. Under a softmax prediction head, we show that accurate token prediction organizes token embeddings according to similarities between the distributions of contexts in which different token types appear, as measured by Hellinger distance, with explicit errors governed by prediction accuracy and token frequency. Meanwhile, the contextual representation provides a low-dimensional coordinate for the conditional distribution of the target token relative to these embeddings. We further introduce a self-consistency principle showing that repeated applications of a shared representation block can progressively refine the contextual representation without introducing additional block parameters. Among representations with the same prediction accuracy, this recurrent construction favors those that can be stably reconstructed from their contexts. Finally, we establish downstream guarantees for token generation, token community recovery, and classification by a linear probe, showing how prediction accuracy and recovered geometry translate into performance beyond the pre-training objective. Together, these results explain how the simple objective of predicting tokens can recover semantic geometry and produce broadly useful representations. A controlled simulation illustrates the theoretical mechanisms.

This paper proposes a novel data-driven algorithm to approximate the dominant eigenfunctions (aka.~modes) of the Koopman operator of nonlinear dynamical systems using neural networks. The relevance of learning the dominant Koopman modes is to approximate nonlinear dynamics by linear ones in a lifted space, thereby enabling simplified control and analysis. To fight the curse of dimensionality arising from using expressive templates (here neural networks) for the mode approximation, the proposed method leverages a power-iteration scheme that directly learns the dominant Koopman modes without explicitly constructing the projection of the Koopman operator on the template of functions. Our approach connects to other approaches in the literature that avoid the curse of dimensionality by learning small dictionaries of functions, but differs from them in that we do not require ``anti-collapse mechanisms'' to ensure that the learned dictionary is expressive enough to approximate the Koopman operator since our power-iteration scheme is designed to converge toward the dominant modes of the projected Koopman operator. The approach is fully data-driven, requiring only sampled state transitions. Theoretical guarantees are provided, showing convergence under increasing sample size and network width (in connection with the neural tangent kernel theorem). Numerical experiments demonstrate that the method achieves accurate and smooth approximations of dominant modes while avoiding the limitations of traditional techniques such as extended dynamic mode decomposition.

Whittle index policies offer a scalable method for restless multi-armed bandits, but under partial observability even determining the indifference subsidy at a single belief requires solving an infinite-horizon belief-state problem with no closed-form value function. Liu [10] addresses this difficulty by linearizing the unknown decision boundary, leading to a linear system and a closed-form approximate Whittle index. However, the resulting threshold uses only a one-step active--passive comparison and does not account for longer-horizon continuation values. We extend this framework to a \emph{$t$-step lookahead threshold policy}. For each subsidy $m$, the threshold is defined by the active-minus-passive advantage under $t$-step finite-horizon value iteration. At $t=1$, the threshold is $m$-independent and recovers the linear threshold of Liu [10]; for $t>1$, it becomes subsidy-dependent through the induced first-crossing structure and tracks the exact decision boundary more closely. The proposed algorithm does not require indexability as an input and includes an indexability verification. Under the original Whittle indexability, we prove that the $t$-step approximate Whittle index converges geometrically to the exact Whittle index, $$|\widehat W_t(ω)-W(ω)|=O(β^t).$$ Numerically, all 2,715 tested three-state instances are verified as indexable according to the proposed criterion. The P95 index error decreases from $2.18\times10^{-2}$ at $t=1$ to $8.93\times10^{-4}$ at $t=8$. In an exact-comparable instance with $β=0.9999$, $t=2$ already recovers the exact Whittle-index ordering. Moderate-depth threshold policies also outperform the one-step baseline and remain close to the optimal dynamic-programming benchmark, while runtime grows mildly with $t$.

Neural PDE training yields a finite checkpoint archive, yet its logged energy errors are inaccessible without the exact solution, while loss-based selection does not necessarily recover the logged energy oracle. For admissible neural approximations of symmetric coercive variational problems, we introduce a reference-free selection rule based on minimizing a computable conforming Riesz monitor. The exact residual-energy identity and conforming projection make the monitor an unconditional lower bound converging monotonically to each logged energy error under nested conforming refinement; under saturation, hierarchical enrichment yields a computable upper estimate and hence a lower-upper bracket. A key finding is that archive selection is order-sensitive: unresolved checkpoint-dependent components can reverse the oracle-non-oracle ranking at finite resolution, so checkpointwise recovery alone is insufficient. For finite archives, we prove uniform recovery, yielding convergence to the logged-oracle error and, without saturation, logged-oracle selection at sufficiently fine auxiliary resolution. Under saturation, the bracket gives a computable near-oracle bound and certifies unique logged-oracle selection upon interval separation. We also bound logging-resolution loss and certify oracle inclusion over prescribed comparison trajectories. The resulting criterion replaces inaccessible exact-error minimization by computable, training-independent post-training selection on the intrinsic energy-error scale, requiring only the computed candidates and the variational problem. Experiments on diffusion and elasticity, including a non-manufactured perforated plate, demonstrate energy-scale calibration, oracle-level selection, and modest post-processing cost.

We develop the Sparse Orthogonal Regression Technique (SORT), a sparse spectral framework for learning orthonormal-basis expansions from noisy and irregularly sampled data. SORT estimates expansion coefficients directly from observations using L1-regularized regression, avoiding explicit quadrature or analytic inner-product evaluation. The central application is data-driven discovery of ordinary differential equations: vector fields are represented in chosen orthogonal bases and learned as sparse coefficient expansions. This provides a complementary route to symbolic regression, grammar-based discovery, and SINDy-style sparse identification by first recovering a compact spectral representation, which can later guide searches for simpler analytic forms. Across the dynamical-system experiments, SORT matches or improves upon library-based sparse-regression baselines when the basis is well adapted to the problem, and shows more stable degradation under sparse sampling, noisy derivative estimates, and representation mismatch. Specific examples illustrate why this representation is useful: if a finite library misses the problem-specific nonlinearity, the resulting model can fail. SORT is not immune to mismatch, but it shifts the problem away from brittle selection among generic terms to basis design adapted to the problem domain. The experiments also show that dominant low-order coefficients persist as model order increases, supporting order-consistent model growth. Beyond equation discovery, the same learned expansion supports nonlinear approximation and estimation of complex, high-dimensional integrals by coefficient readout. Overall, SORT provides a reusable intermediate representation for system identification, approximation, and integration, while making basis design an explicit part of the scientific modeling problem.

Learning mappings between infinite-dimensional objects is a central challenge in scientific machine learning. We introduce a general kernel-based encoder-decoder framework for operator learning that separates observation, representation, learning, and reconstruction. We develop this framework for multi-input, multi-output operator learning, where operators map between products of potentially distinct function spaces. Our approximation theory shows that, although the number of inputs and outputs can increase, the convergence rate is governed by the most challenging constituent approximation problem rather than the overall problem dimension. The framework leads to practical kernel methods with closed-form training and inference, combining mathematical tractability with computational efficiency. We further specialize the approach to multiple operator learning by introducing KernelMO, a family of kernel methods with complementary operator-valued and product-space formulations. Across five families of parametric partial differential equations, the proposed methods achieve competitive or state-of-the-art predictive accuracy while reducing training and inference costs relative to neural operator architectures and deep learning based models, offering an efficient and lightweight alternative.

Many reinforcement learning systems, from fleet management to traffic signal control, must serve an objective that changes dynamically after deployment, and retraining a policy for each new objective is prohibitively expensive. For a single agent, this problem is well understood: successor features with generalized policy improvement, together with their universal extension, recombine a library of learned policies into a policy for any new objective, with a guarantee that the result is never worse than any policy in the library. However, multi-agent transfer has received far less attention, and the common practice of letting each agent recombine its own library independently inherits the recipe but not the guarantee. We prove that this independent composition can produce joint behavior strictly worse than every policy in the library, because recombining teammates changes the environment each agent faces and invalidates the values it relies on, a failure with no single-agent counterpart. We further show that the only unconditionally safe fixed rule is synchronized composition, which moves the whole team to one jointly trained policy but cannot serve objectives that assign different goals to different agents. To attain safety and flexibility at once, we propose MA-USFA, a hierarchical method with two layers: a lower layer of universal successor feature approximators that predicts each agent's successor features while conditioned on its teammates' objectives, and an upper composer that selects, across agents, which library entry each agent should follow and supplies the cross-agent correction a per-agent value cannot represent. Trained once over the distribution of objectives, it is applied at deployment with no per-task adaptation.

Clustering is a fundamental class of data analysis techniques with the most important representatives being centroid-based methods like $k$-means. Such methods are strongly connected to quantization problems, which aim to approximate general probability measures with discrete ones. For example, $k$-means corresponds to quantization with respect to the Wasserstein distance. While Wasserstein quantization clusters points within a fixed space, this paper studies Gromov-Wasserstein (GW) quantization, which additionally aims at clustering the ambient geometry of the space. We show existence of solutions to the GW quantization problem and give a characterization that justifies an analogue to the $k$-means algorithm (Lloyd's algorithm) to approximate them numerically. We further calculate the quantization rate for usual Euclidean geometries that are used in the GW context, and relate it to standard Wasserstein quantization rates. Finally, numerical experiments show that GW quantization opens up many modeling possibilities beyond normal clustering methods (e.g., for geodesic distances of 3D shapes or structured pruning of neural networks) and that the introduced algorithm leads to useful numerical solutions with approximation quality often in line with theoretically optimal rates.

We study the approximation of dynamical systems by semi-autonomous neural ordinary differential equations (SA-NODEs) over long time horizons. For a single network trained on the whole horizon, the available error bound deteriorates double exponentially in the horizon length. We develop two training strategies that avoid this barrier, each built on a reset of the state. The model predictive strategy partitions the horizon adaptively and restarts every window from observed data: when training meets a prescribed tolerance on every window, the composite model meets it uniformly in time, with a parameter budget linear in the horizon for targets with a bounded, uniformly regular reachable tube. The Floquet strategy addresses autonomous targets with a stable limit cycle and uses no data at deployment: a certified contraction of the learned return map confines the error to linear growth in the number of elapsed periods. For the time-periodic architecture we deploy, the scalar certificate degenerates; we prove instead a uniform-in-time orbital guarantee whose hypotheses are measured on the trained model, and an obstruction showing that, for an exactly periodic learned field, small one-period error and a contracting stroboscopic map cannot hold at once. Numerical experiments on four benchmarks confirm the predicted error laws and measure the hypotheses of every guarantee.

Reducing the number of focal elements of a mass function is classically driven by an intrinsic distance, such as Jaccard or Jousselme, that keeps the approximation close to the original as a body of evidence. We consider instead the case where the mass function feeds a linear combinatorial optimisation problem with evidential costs. What should then be preserved is not the closeness of the two mass functions, but the quality of the decision they induce. We introduce a decision-aware approximation that targets the regret of the decision: one decides with the cheaper approximation and is evaluated under the true mass function. On a minimal shortest path, the distance-optimal approximation flips the decision while a decision-aware merge preserves it, and this occurs on a non-negligible fraction of random instances. We prove a one-point bound that localises the regret at the true optimum, turn it into an exact dynamic program for the scalar case, and extend it to an online version that prunes focal elements before the final cost is known. In experiments the decision-aware compressor flips the decision less often than representation-aware compression, for both the linear criterion and a non-linear proxy read-out.

Universal approximation in reservoir computing is typically associated with a class of reservoirs. We show that universality can be associated with a single reservoir, considering a minimal setup of a single atom in front of a mirror. In its linear-transducer limit, our reservoir is a universal approximator of fading-memory maps under an operating class of checkable conditions, with a rate constant measured at the operating point. A given reservoir can reach arbitrary accuracy by changing measurement settings. The proof gives an explicit recipe: for a target accuracy, it specifies the required physical resources and resonator modes. Enlarging the number of accessible modes increases the matchable kernel span without reducing capability. Beyond the linear limit, the atom's saturation replaces high-order polynomial readouts, and the device operates on real-world tasks alongside classical baselines. Our results highlight an example of universality with a minimal quantum setup.

How expressive is prompting a transformer? Answering this question is important for separating the roles of prompting, architecture, and pretraining in transformer models, and for determining whether task-specific behavior must be stored in model weights or can instead be induced at inference time through the prompt. We show, in an approximation-theoretic sense, that pretraining is optional: a single-layer softmax attention network with random, untrained weights can approximate any Hölder function on a compact manifold when steered by an appropriate soft prompt. Guided by the connection between softmax attention and kernel methods, we construct explicit soft prompts (a prompt per target function, independent of the query) as solutions to linear systems matching attention logits to Gaussian kernel exponents, under which the frozen transformer emulates the classical Nadaraya-Watson kernel estimator. The construction requires only a mild rank condition on the weights, which we show holds almost surely under Gaussian initialization. The prompted network inherits the theoretical guarantees of kernel regression, leading to universal approximation theorems with minimax-optimal rates that depend on the intrinsic dimension. We further quantify the cost of prompting, exposing a tradeoff between the norm of the constructed soft prompt tokens, prompt length, and hidden dimension. Numerical experiments corroborate the constructions and predicted rates.

In this paper we provide Monte Carlo and deep neural network approximations for stochastic representations of solutions to linear elliptic partial differential equations with constant diffusion, drift and killing. Building on the modified Walk-on-Spheres algorithm of Beznea et al. (arXiv:2209.01432), we introduce Monte Carlo estimators that explicitly incorporate sampled random times arising in the analyzed stochastic representations. We establish uniform error bounds for these estimators and show that, under suitable assumptions, a prescribed approximation accuracy is achieved with sample complexities growing at most polynomially in both the inverse accuracy and the problem dimension. Furthermore, we prove a deep neural network approximation result for the stochastic representations. Assuming suitable neural network representations of the boundary data and the distance function to the boundary, we use the constructed Monte Carlo to design deep neural networks that approximate the representation uniformly with a number of parameters growing at most polynomially in the inverse accuracy and the problem dimension. These results extend previous complexity analyses to a broader class of elliptic equations involving drift and killing.

In this paper we develop quantitative approximation results for shallow neural networks constructed using a dictionary based on metaplectic operators. First, we extend the concept of Barron spaces by considering a symplectically motivated extension of the Fourier transform, known as the metaplectic transform. Then, after establishing embedding between metaplectic Barron spaces and Sobolev spaces we consider a neural metaplectic dictionary and we prove Monte-Carlo approximation bounds for metaplectic Barron functions using finite linear combinations of atoms of the dictionary. Finally, we validate the introduction of the neural metaplectic dictionary by devising a deep neural network architecture that uses as building blocks the atoms of the dictionary. We test it to approximate solutions of time-dependent Schrödinger equations, demonstrating better performance compared to classical phyisics informed neural networks architectures.

We introduce conditional cylindrical neural networks for approximating functionals of conditional laws in McKean-Vlasov equations with common noise. Fourier moments of the initial law and truncated signatures of the time augmented common noise are mapped by a mixture density network to a Gaussian mixture approximation of the conditional law. A cylindrical neural network then evaluates the target functional through analytic integrals against this predicted measure. Rough path well posedness and stability provide a conditional law map that is continuous in the initial distribution and the rough driver and agrees almost surely with the classical conditional law at the Itô Brownian lift. Combining this continuity with Fourier separation, signature uniqueness, Wasserstein density of Gaussian mixtures, and neural universal approximation, we prove an $L^2$ universal approximation theorem for continuous square integrable functionals. The numerical study implements the resulting two stage procedure on six examples, including non Gaussian initial laws, nonlinear drift, multiplicative common noise, and a two dimensional state. Independent particle references are used when no closed form law is available. The learned conditional law and functional approximations consistently improve on the empirical particle plug in, and additional experiments examine feature sensitivity, training from one terminal observation per common noise scenario, and Itô--Stratonovich consistency.

We develop a rigorous theory of discrete residual least-squares approximation for elliptic spectral equations $\mathfrak L_βu=f$ using linearized ReLU$^k$ neural networks on the sphere, where $\mathfrak L_β$ is a positive elliptic spectral multiplier of order $β$. Given a parameter set $Θ_n=\{θ_{j}^*\}_{j=1}^n\subset\mathbb S^d$, we approximate $u$ in the linearized network space $L_n^k(Θ_n)$ by the discrete residual on the collocation points $\{η_i^*\}_{i=1}^m$ \begin{equation*} u_{n,m}\in\arg\min_{v_n\in L_n^k(Θ_n)}\frac1m\sum_{i=1}^m\left(f(η_i^)-\mathfrak L_βv_n(η_i^)\right)^2. \end{equation*} With $k>\frac{d-1}{2}+β$, for antipodally quasi-uniform network parameter sets and any quasi-uniform collocation points with $m\gtrsim n$, we prove that \begin{equation*} |u-u_{n,m}|{\mathcal H^β(\mathbb S^d)}\eqsim|f-\mathfrak L_βu{n,m}|{\mathcal L^2(\mathbb S^d)}\lesssim n^{-\frac{r}{d}} \begin{cases} |f|{\mathcal W^{r,p}(\mathbb S^d)},&\frac{d}{p}<r\leq \frac{d}{2},~p>2,\ |f|{\mathcal H^r(\mathbb S^d)},&r>\frac{d}{2}. \end{cases} \end{equation*} We also establish a high-probability residual estimate, up to a logarithmic factor and an arbitrarily small smoothness loss, for i.i.d.\ uniformly distributed collocation points. The key analytical ingredient is a Bernstein inequality for linearized ReLU$^k$ network spaces. If $\underline h$ denotes the antipodal separation distance of the network parameters, then \begin{equation*} |v_n|{\mathcal H^r(\mathbb S^d)}\lesssim\underline h^{-(r-s)}|v_n|_{\mathcal H^s(\mathbb S^d)},\qquad 0\leq s<r<k+\tfrac12. \end{equation*}

Deep Operator Networks (DeepONets; arXiv:1910.03193) typically encode an input function through point values on a fixed discretization. Building on the Topological DeepONet framework of Ismailov (arXiv:2603.11972), we replace point samples by continuous linear functionals drawn from the continuous dual of a Hausdorff locally convex space $({V},\{p_α\}_{α\in A})$, whose topology is generated by a point-separating family of seminorms rather than a single norm, and develop fixed and adaptive functional measurement systems. Measurements are combined with the coefficient-space Two-Step procedure of Lee and Shin (arXiv:2309.01020), while a training-only decoder and regularization stabilize the adaptive coordinates. We derive a discrete error decomposition separating measurement, output-basis, and neural-approximation errors, together with a Barron-rate refinement. The framework is evaluated on the antiderivative operator, a non-normable locally convex input space, heterogeneous Darcy flow, a controlled operator, and fixed-time and time-evolving Navier-Stokes vorticity operators. In the heterogeneous Darcy problem, the functional models retain nearly resolution-independent errors of 5.5-5.6% on unseen grids, while in the controlled problem adaptive measurements reduce the mean error below 1.2%. For the fixed-time Navier-Stokes problem, the Adaptive Topological DeepONet is the most accurate DeepONet-based model, attaining a mean relative $L^2$ error of 1.685% +/- 0.017% using 128 functional coordinates. A comparably sized Fourier neural operator (FNO; arXiv:2010.08895) achieves the lower error 0.832% +/- 0.172%, but requires the full 64x64 input field, twice the training time, and 10.7x greater peak GPU memory. The formulation provides compact, interpretable, and discretization-portable coordinates in the continuous dual $V'$, including for non-normable input spaces.

We present a unified theoretical framework connecting classical universal approximation theory, Fourier-feature approximation, and diffractive optical processors. We show that phase-encoded diffractive processors implement finite Fourier-feature expansions whose mathematical completeness follows from Fourier/Stone-Weierstrass arguments, while their physical realizability is governed by finite coefficient synthesis through optimized spatially varying coherent point-spread functions (PSFs). Our analyses derive approximation-error bounds that separate Fourier truncation, PSF-synthesis, input phase error, optical hardware, readout, and noise contributions; establish scaling relationships linking approximation complexity to optical degrees of freedom and input/output space-bandwidth products; derive photon-budget and throughput limits imposed by photon statistics; formulate finite-class statistical learnability bounds for phase-quantized diffractive function approximators; and analyze the impact of spatially incoherent illumination. We further analyze coherent optical cascadability and show that quadratic feature expansion through coherent mixing and optical readout provides a mechanism for enhanced representation while remaining fundamentally distinct from the depth-separation results established for digital neural networks. Our analyses provide a rigorous theoretical foundation for diffractive nonlinear function approximation and establish quantitative relationships among mathematical expressivity, optical hardware resources, statistical learning, and physical performance limits, thereby offering general design principles for large-scale analog optical computing systems.

Robotic planning often involves multiple objectives with complex priority relationships, such as safety, efficiency, and regulatory compliance. Rulebooks formalize these relationships, allowing partial ordering of objectives that generalizes both Pareto and lexicographic dominance. Computing the full set of rulebook-optimal solutions, however, is computationally expensive. To address this challenge, we introduce the concept of epsilon-rule-dominance, a principled notion of approximate dominance under rulebooks, and propose RApex, a best-first search algorithm that efficiently computes a compact set of epsilon-approximate rulebook-optimal solutions. RApex leverages dimensionality reduction, a technique used to speed up existing multi-objective search algorithms, while respecting rule hierarchies by maintaining separate closed sets and performing dominance checks over truncated and residual rule sets. We provide a formal analysis of RA*pex, proving that every rulebook-optimal solution is epsilon-rule-dominated (a generalization of approximate dominance we introduce) by at least one solution in the returned set. Empirical results demonstrate that our approach achieves computation times over two orders of magnitude faster than existing methods.

We propose a family of random feature maps for scalable kernel machines on low-dimensional subspaces, ie on the Grassmannian manifold. Such representations are useful when data classes or clusters are well described by the span of a few samples. Classical Grassmannian kernels, including the projection and Binet-Cauchy kernels, require full Gram matrices, which leads to prohibitive computational and memory costs for large high-dimensional subspace datasets. We address this limitation using random features based on rank-one projections of subspace projection matrices followed by bounded non-linear transforms, either periodic or binary, to control the resulting distributions. We show that inner products in the random feature space approximate well-defined rotation-invariant Grassmannian kernels that depend only on the principal angles between subspaces. When the number of features is sufficiently large relative to the intrinsic subspace dimension, the approximation holds uniformly over all fixed-dimensional subspaces with high probability. For periodic transforms, the approximated kernel has a closed-form expression with tunable behaviour between inverse Binet-Cauchy and Gaussian-type regimes. Binary transforms yield compact one-bit subspace features, although no closed-form kernel is known. Structured rank-one projections based on randomised fast Fourier transforms further reduce computation without sacrificing practical accuracy. Experiments on synthetic data and ETH-80 classification tasks show that these features accurately preserve Grassmannian geometry while reducing computation, memory, and storage. Rank-one embeddings therefore provide a practical and scalable alternative to classical Grassmannian kernels.

Deep reinforcement learning (DRL) approaches for flexible job shop scheduling (FJSP) heavily rely on attention-centric architectures to achieve state-of-the-art performance. However, these models suffer from excessive parameter counts and prohibitive inference latency as problem scales expand. While liquid neural networks (LNNs) offer a parameter-efficient alternative for modeling adaptive state evolution, their inherently sequential dynamics bottleneck computational efficiency. To resolve this trade-off, we propose PLAN (Parallel Liquid-inspired Approximation Network), a lightweight representation learning framework that reformulates continuous liquid-state dynamics into a discretized and parallelizable formulation. PLAN structurally decouples state evolution from context aggregation, where liquid-inspired updates handle the primary evolving state representation, and a lightweight context aggregation module provides complementary global context. Furthermore, PLAN acts as a versatile, plug-and-play backbone that generalizes to complex FJSP variants, pairing with a compact stochastic module for stochastic FJSP and replacing heavy heterogeneous graph transformers in multi-faceted dynamic FJSP. Extensive evaluations across deterministic, stochastic, and multi-faceted dynamic FJSP benchmarks show that PLAN reduces the average makespan by 1.2%, 1.4%, and 2.3%, respectively, compared with the corresponding state-of-the-art baselines, with the improvement reaching 10.2% in one benchmark setting. PLAN also reduces average inference latency by 13.2%, 31.7%, and 26.9%, respectively, with a maximum reduction of 69.2% on the largest instances, while using only 22$-$47% of the baseline parameters.

We present a simple Gaussian approximation to the finite-sample distribution of the classical ridge regression estimator. Our approximation captures the fact that, in finite samples, the ridge regression estimator trades off bias and variance to reduce estimation and prediction error. Our approximation is based on nonstandard asymptotics where $i)$ we let the estimator's regularization parameter grow proportionally to the sample size; and $ii)$ we treat the population regression coefficients as \emph{local} to the reference vector that defines the estimator's direction of shrinkage. In contrast to other asymptotic approximations in the literature, we allow for general forms of heteroskedasticity and autocorrelation in the data generating process (at the cost of considering a low-dimensional model where the number of covariates is not allowed to grow with the sample size). We use our simple Gaussian approximation to propose two new strategies to select the regularization parameter for the ridge regression estimator. The suggested strategies select the regularization parameter to minimize either average or worst-case excess prediction risk, where risk is computed using our suggested Gaussian approximation.

In this paper, we consider the parameter synthesis and optimization problem for parametric Markov decision processes (pMDPs), the extension of classical MDPs where exact probability values are replaced by parametric expressions. Computing the rational function $f_{\lsf}$ that maps parameter valuations to the satisfaction value of a PRCTL property $\lsf$ is a computationally expensive task, particularly for pMDPs where the optimal policy may vary across the parameter space. We adopt the \emph{scenario approach} to efficiently synthesize a probably approximately correct (PAC) approximation $\ApproxFunOfProperty{f}$ of $f_{\lsf}$: by sampling parameter configurations and solving a linear program, we obtain a polynomial approximation whose error margin $\margin$ is guaranteed, with prescribed confidence, for all but an $\errorRate$-fraction of the parameter domain under the sampling distribution. We further show how this PAC framework can be combined with statistical model checking (SMC), enabling the analysis of black-box parametric models. Building on the PAC approximation, we integrate the DIRECT (DIviding RECTangles) algorithm for derivative-free global optimization over the parameter space. We establish conditional optimality-gap guarantees: under explicit Lipschitz and PAC-good-set assumptions, the difference between the true optimum $f_{\lsf}(\parameters^{*})$ and the value found by DIRECT is bounded by a partition-diameter term and, in the PAC case, an additional approximation-error term. An empirical evaluation on 2997 benchmarks focuses on the new DIRECT-based optimization component. The results show that DIRECT variants solve fewer instances than the scenario optimizer, but on their common successful instances they often return slightly better objective values and usually run faster, while remaining close to the scenario values within the PAC margin.

Traditional approximation theory measures convergence rates in terms of the number of parameters or degrees of freedom. However, practical computation operates under finite precision: parameters must be encoded using a finite number of bits. Therefore, approximation efficiency should be evaluated in terms of computational bit complexity, which is intrinsically connected to the metric entropy of the underlying function class. In this work, we develop a unified approximation framework based on binary encoding and metric entropy. We analyze classical methods (including polynomial approximation, sparse grids, and finite elements) as well as shallow and deep neural networks, and compare their approximation rates for function classes with comparable metric entropy. We observe that, when evaluated in terms of bits, most classical methods are in general suboptimal relative to the intrinsic limits dictated by metric entropy, while neural network methods may exhibit different behaviors. We show that when complexity is measured in bits rather than parameters, no method fundamentally exceeds the approximation order achieved by classical approaches. Our results also indicate that many seeming advantages of neural networks, including dimension-independent rates and superconvergence phenomena, stem from differences in function class complexity rather than intrinsic architectural superiority. In this sense, the traditional curse of dimensionality can be misleading; the fundamental limitation is instead a curse of bit complexity, governed by metric entropy.

We establish a dimension-efficient neural network approximation theory for solutions to fractional parabolic equations with lower-order drift and potential terms. By introducing anisotropic spectral Barron spaces, which measure temporal and spatial regularity separately in frequency space, we first develop a dimension-independent maximal regularity theory for these equations, using dimension-independent multiplication estimates and the method of continuity to incorporate the lower-order terms. A key technical novelty is the application of the Vandermonde matrix to the global-in-time extension of the finite-time fractional heat semigroup with sufficient regularity at the initial time, thereby enabling analysis of the forward-in-time evolution via the global space-time Fourier structure of anisotropic Barron norms. We also show that a corresponding uniform-in-time estimate of the spectral Barron regularity generally fails. Finally, we derive $n^{-1/2}$ two-layer approximation bounds in mixed Sobolev norms for non-constant periodic activations and, under additional anisotropic Barron regularity, for non-periodic activations satisfying a polynomial-decay condition.

Engression (Shen and Meinshausen, 2024) learns a conditional distribution by fitting a generative model $Y = f(X,\varepsilon)$ under the energy score, a strictly proper scoring rule. We provide a theoretical error analysis of engression implemented with deep neural networks. We decompose the excess risk into three components: the approximation error, the stochastic error, and the Monte Carlo error. Based on this decomposition, we establish convergence rates under the assumption that the target conditional generator admits a compositional smoothness structure.

Accurate option prices do not imply accurate recovery of the latent risk-neutral density. We study this distinction with two complementary benchmarks. A controlled benchmark exposes simulator-truth densities for latent evaluation, while a chronological NIFTY benchmark tests only held-out market prices. A two-component lognormal mixture has the lowest aggregate price, $L^1$, Wasserstein, and fixed-tail errors on the synthetic benchmark. Learned operators retain narrower strengths: DeepONet reduces 1% quantile and variance error by 39.0% and 34.6% relative to the mixture, and a quote transformer reduces $L^1$ by 16.4% on the structurally misspecified Merton family. A numerical conditioning analysis explains why these rankings can differ: after enforcing mass and forward constraints, 95 of 126 pricing directions are numerically null, and two densities separated by $L^1 = 0.061$ produce identical prices on the covered strikes. On 524 held-out NIFTY calls, validation-selected test-time adaptation reduces DeepONet RMSE by 28.3%, but per-expiry mixture and SVI fits remain much more accurate. The evidence supports target-dependent inductive bias, not a universal winner.

We investigate message-passing graph neural networks with random node features. Random node features are known to enhance the expressiveness of graph neural networks (GNNs) both theoretically and empirically. Here, we establish a novel universality result focusing on permutation-equivariant neural networks (PENNs), a class of GNNs built from feedforward neural network components that subsumes many prominent GNN architectures. We show that PENNs, combined with partially random node features, can approximate arbitrarily well in probability any measurable permutation-invariant or permutation-equivariant function on directed graphs of fixed size with multidimensional node and edge features. For $k$-times continuously differentiable functions, $k\geq 2$, we also derive upper bounds on the approximation rates, relating the complexity of the feedforward components of a PENN in terms of layer depth and number of nonzero weights to the desired approximation accuracy.

The $k$-means++ algorithm is a standard and widely used seeding method for $k$-means clustering, but for a fixed number $k$ of centers its worst-case expected approximation ratio is $Θ(\log k)$. We consider the same algorithm when an adversary first fixes the dataset and some $K$; the number of centers $k$ is then chosen uniformly from $\{K,\ldots,2K-1\}$. We prove that $k$-means++ is an $O(1)$-approximation with constant probability in this budget-smoothed setup.

Matricized low-rank approximation via SVD is a standard surrogate for tensor decompositions, but entry-wise reconstruction error fails to capture multiway geometric degradation. Under an orthogonal Tucker model, we characterize this degradation using two metrics: cross-mode Direction Loss, measuring geometric subspace deviation from rank truncation and noise rotation, and Interaction Loss, quantifying multilinear interaction distortion in the core tensor. We prove that squared relative reconstruction error orthogonally decomposes into interaction loss and out-of-subspace energy, and derive a Wedin-type bound establishing the stability of a plug-in Direction Loss estimator. Experiments on synthetic and hyperspectral datasets demonstrate that nearly identical reconstruction errors can yield markedly different structural-loss profiles; hyperspectral patches with comparable reconstruction errors exhibit up to a 4.6-fold difference in Direction Loss, correlating with severe visual blurring.

Fully homomorphic encryption (FHE) provides strong cryptographic guarantees for private inference, but deploying transformer models under FHE remains prohibitively expensive. A key bottleneck is that non-linear operations such as softmax, normalization, and activation must be replaced with polynomial approximations compatible with the CKKS scheme, and the multiplicative depth consumed by these approximations dominates inference cost. Recent frameworks have advanced approximation techniques, yet all rely on manually configured approximation hyperparameters (e.g., number of iterations, polynomial degree), applied uniformly across all layers. While convenient, this uniform-configuration approach is overly rigid: different layers can tolerate different levels of approximation error without degrading predictive accuracy, and uniform configurations cannot exploit this variability to reduce latency. Allowing each layer to adopt its own configuration, however, causes the search space to explode with model depth, reaching roughly $10^{84}$ configurations for BERT/ViT (12 layers) and $10^{225}$ for LLaMA3 (32 layers), rendering manual exploration practically impossible. We present ATLAS, an automated framework that configures per-layer approximation settings by formulating the problem as a multi-objective optimization over latency and predictive accuracy. The resulting problem is inherently difficult: 1) competing objectives over a large decision space (120 or 320 variables for BERT/ViT or LLaMA3); 2) expensive evaluation, as each configuration takes 70-1,000 seconds even in cleartext; and 3) sparse optimization signals, as 35-50% of candidate configurations yield numerically invalid solutions. ATLAS addresses these challenges through a two-stage optimization strategy that progressively relaxes layer-wise constraints, combined with surrogate models to accelerate evaluation.

In this paper, we study the extension of Neural Jump ODEs to infinite-dimensional function spaces. In particular, the underlying process $X$ now takes values in $L^2(Ξ, \mathbb{R}^{d_X})$ instead of $\mathbb{R}^{d_X}$ and the Operator NJ-ODE approximates the optimal predictor of this process by producing a representative of the conditional expectation. The NJ-ODE model is a framework for online learning the optimal prediction of continuous-time stochastic processes, given discrete, possibly irregular and incomplete past observations. In a series of works, this model has been extended to deal with generic path-dependent processes, with observation noise and dependent observations, with long-term predictions, and with input-output systems. However, throughout all of these works, the underlying processes were restricted to be finite-dimensional. In particular, function-valued problems, like yield curve or volatility surface predictions, could only be handled through discretization, which inherently leads to a loss of information. In this work, we build on ideas from Neural Operator methods that allow us to extend the NJ-ODE framework to an infinite-dimensional output process. To prove convergence of the NJ-ODE to the optimal prediction process, we develop a new approximation strategy that also generalizes previous works in the finite-dimensional setting by considerably weakening the assumptions.

Score-based diffusion models have achieved remarkable empirical success in generative modeling, yet their approximation-theoretic foundations remain incomplete. In particular, although classical universal approximation theorems guarantee that neural networks can approximate score functions, it remains unclear whether such approximation guarantees translate into approximation of the probability distributions generated by reverse diffusion processes. In this paper, we establish a rigorous quantitative connection between these two notions. Specifically, we prove that if a neural network approximates the true score function sufficiently accurately, then the probability distribution generated by the corresponding reverse diffusion model is close to the target data distribution in Kullback-Leibler (KL) divergence, up to an irreducible mismatch between the terminal distribution of the forward diffusion process and the prior used to initialize the reverse process. More precisely, we derive an explicit upper bound on the distribution approximation error in terms of the score approximation error, the diffusion noise schedule, and the terminal prior mismatch. Our analysis combines Hornik's universal approximation theorem, Girsanov's theorem on path space, and the data processing inequality for relative entropy. Complementary to recent work that studies score approximation under finite-sample statistical settings and structural assumptions on the data distribution, our work develops an approximation-theoretic analysis based on classical neural network approximation theory. The resulting theorem provides a simple and explicit guarantee linking neural network approximation of score functions to approximation of the probability distributions generated by reverse diffusion models.

Motivated by the numerical computation of the Mean Escape Time (MET) $τ:Ω\to\mathbb{R}$ of a stochastic process from a bounded domain $Ω\subseteq\mathbb{R}^d$, we study elliptic Dirichlet boundary value problems (BVPs) using boundary-enforced Physics-Informed Neural Networks (PINNs), in which the Dirichlet condition is imposed exactly by multiplying the network output with a predefined distance-to-boundary approximation $ρ$. Combining approximation-theoretic and statistical-learning arguments for Rectified Quadratic Unit (ReQU) and hyperbolic tangent (tanh) networks, we derive a priori error bounds that make explicit the dependence on $ρ$. In particular, we show that exact boundary enforcement alone is not enough for $H^2(Ω)$ error bounds, and that a sufficient and essentially necessary condition is for $ρ$ to be a smooth distance approximation $\textit{normalized to first order}$, of the kind constructed in arXiv:2104.08426 [math.NA]. We thereby identify this subclass of $\textit{boundary-adapted}$ PINNs as the appropriate neural network ansatz for solving Dirichlet BVPs. Numerical experiments support the theory, showing that appropriate choices of $ρ$ improve accuracy and convergence, while poorly chosen distance functions can substantially degrade the solution. Our proof also yields new VC-dimension bounds for hypothesis spaces of higher-order derivatives of ReQU and tanh networks, together with new approximation bounds for shallow ReQU networks in higher-order Sobolev norms, all of which are of important independent interest.

The Universal Approximation Theorem states that a neural network with a single hidden layer is sufficient to approximate any continuous univariate function on a compact domain to arbitrary error. However, the uniqueness of such neural network representations is not guaranteed, raising questions about practical identifiability. In this work, we address this concern by analyzing functional equivalence and geometric diversity of neural network approximations to a few elementary mathematical functions. The analysis includes an extensive study of single-layer neural networks and multilayer perceptrons under noisy and noise-free conditions. Beyond just network capacity, we study the geometric properties through the lens of sloppiness, characterized by the eigen spectrum of the Hessian of the cost function and the effective rank to quantify the dimensionality of parameter space. The study reveals large equivalence classes of functionally indistinguishable yet geometrically diverse networks that consistently exhibit low effective rank and structural redundancy. Finally, a model select criterion is proposed for identifying optimal models based on parsimony, ease of estimation, and inference efficiency.

The paper presents a new parameter-efficient adaptation method called ChebyMA (Chebyshev Manifold Adaptation). ChebyMA adopts weight matrices through a multi-surface superposition of Chebyshev polynomial bases evaluated on learnable coordinates and combined via trainable coefficient matrices, replacing standard linear projections with highly expressive continuous function approximation. Theoretically, we establish an Approximation Expressivity Theorem, proving from the perspective of function approximation theory that single-manifold ChebyMA guarantees convergence in Frobenius norm error of reconstruction. Besides, drawing on Kolmogorov $n$-width intuition, we demonstrate the expressive advantages of multi-manifold superposition ($S > 1$) in decoupling high-dimensional complex features. Experimental results on Computer Vision CIFAR datasets(CIFAR-10, CIFAR-100)\cite{CIFAR} and Natural Language Processing (AG News, SST-2) datasets demonstrate that ChebyMA consistently achieves a superior parameter-accuracy Pareto front compared to standard full-parameter fine-tuning, LoRA\cite{LoRA}, TLoRA\cite{TLoRA}, and StelLA\cite{StelLA}. ChebyMA significantly outperforms other tested methods in tested datasets, validating its solid theoretical foundation for generality with purely vectorized computations.

We develop a non-asymptotic approximation, sampling, and finite-iteration optimization theory for variational physics-informed approximation of uniformly monotone nonlinear multiscale elliptic equations. For boundary-compatible neural feature classes, the population error splits into approximation, empirical quadrature, and projected-gradient terms, with all non-approximation constants uniform in the microscopic scale $\varepsilon$. Assuming a quantitative corrected $H^1$-estimate, a two-scale state class yields $$\mathcal A_m^\varepsilon \le C\bigl(\varepsilon+Φ_{0,m_0}^2+Φ_{1,m_1}^2\bigr)$$ in arbitrary dimension. We further introduce a convex primal-dual physics loss whose population value is a computable upper certificate for the state error. With additional flux-corrector regularity, a divergence-compatible two-scale flux class gives a certified state-flux bound combining $O(\varepsilon)$ approximation, state and flux feature errors, empirical sampling error, and an $O(K^{-1})$ optimization term. In contrast, for general periodic nonlinear fluxes satisfying a natural nondegeneracy condition, the empirical Rademacher complexities of strong-residual and squared-residual classes are bounded below by constant multiples of $(\varepsilon\sqrt N)^{-1}$ and $(\varepsilon^2\sqrt N)^{-1}$, respectively. These optimizer-independent lower bounds hold in every spatial dimension. Numerical experiments confirm the predicted $\varepsilon$- and $N$-scalings for nonlinear fluxes in $d=1,2,3$, validate every computed primal-dual certificate, and show that corrector-enriched classes substantially reduce energy and $H^1$ errors as the microscopic scale is refined.

We address fundamental challenges in representing and computing $\mathbb{R}^{d}$-valued predictable square-integrable processes over $[0,T]$, collected in the space $\mathcal{H}^2_T(\mathbb{R}^{d})$. These processes are central to continuous-time stochastic control, reinforcement learning, and mathematical finance. Although Wiener-chaos expansions offer strong theoretical tools, traditional computational methods are hindered by the need for large chaos dictionaries and high-order iterated integrals. To overcome these obstacles, we introduce NeuralChaos -- a neural operator architecture that produces elements of $\mathcal{H}^2_T(\mathbb{R}^{d})$ using only finitely many evaluations of the driving Brownian motion, while preserving predictability and square-integrability. We prove that NeuralChaos is dense in $\mathcal{H}^2_T(\mathbb{R}^{d})$ and achieves the best $N$-term chaoslet approximation rates for compressible and Malliavin--Sobolev regular processes. Moreover, compressibility is shown to be typical for processes from $\mathcal{H}^2_T(\mathbb{R}^{d})$ under non-degenerate sub-Gaussian sampling. In contrast, we show that finite-dimensional Markovian neural SDE models constitute a meagre and Gaussian-null subset in $\mathcal{H}^2_T(\mathbb{R}^{d})$, regardless of discretization, whereas compressible processes are generic. Numerical experiments on a stochastic optimal control problem and dynamic hedging highlight the practical effectiveness of our approach. Our results enable more efficient and expressive modelling in stochastic analysis and mathematical finance.

We study the contextual combinatorial semi-bandit (CCSB) problem with general reward function approximation. At each round, the learner observes a context, selects a combinatorial action consisting of a subset of basic arms, and receives the reward of each selected arm; the goal is to maximize the cumulative reward over time. We propose SquareCB.Comb, a computationally efficient algorithm that, at each round, solves a convex optimization problem to sample a combinatorial action that balances exploration and exploitation. SquareCB.Comb scales to large arm sets and imposes no structural assumptions on the action set beyond a cardinality bound of $m$ on each combinatorial action. We prove that SquareCB.Comb achieves a minimax optimal regret bound of $O(\sqrt{m A T \log |\mathcal{F}|})$, where $A$ is the number of arms, $m$ is the maximum number of arms in a combinatorial action, $T$ is the time horizon, and $\mathcal{F}$ is the reward function class. In the realizable setting, this bound matches the state-of-the-art regret guarantees achieved by policy search-based algorithms in the more restricted slate recommendation settings, while simultaneously generalizing to arbitrary combinatorial action structures and general reward function approximation.

We develop a convergent scheme to train neural networks involving analytic activation functions based on gradient flows. Convergence properties are guaranteed by Lojasiewicz theory. The main advantage of this approach is its simplicity of implementation. The coefficients of the network are approximated by solving a system of ordinary differential equations. We test the method by constructing residual neural network approximations of solutions of parametric problems. The dependence of the solutions of simple ordinary differential equations on a few parameters is correctly reproduced. The solutions of inverse problems involving wave constraints which depend on a few parameters can be reasonably approximated, even in regions in which the problem is severely ill posed.