NEOPIC: A Neural Operator Framework for Particle-based Kinetic Plasma Simulations
Sriramkrishnan Muralikrishnan, Chelsea Maria John, Jeffrey Kelling, Michael Bussmann, Robert Speck, Paul Gibbon
AI 速览：提出NEOPIC，一种直接作用于粒子经验测度的神经算子，作为即插即用的场求解器，用于粒子动理学等离子体模拟。
We introduce NEOPIC, a framework for learning the field operator in particle-based methods applied to kinetic plasma simulations. Its key novelty is a neural operator that acts directly on the empirical measure of the particles, rather than on a grid, and is evaluated at the particle positions. No mesh or tree is needed and, unlike grid-based neural operators, no deposition onto a mesh takes place, so the particle noise that dominates the inputs of particle simulations is not aliased. Since it replaces only the field solve, it serves as a plug-and-play field solver for particle-based methods: it can be trained on data from particle-in-cell, particle-in-Fourier or tree codes, and combined with any time integrator, including implicit schemes. We prove that the discrete evaluation coincides with the continuous operator at any number of particles (representation equivalence) and derive an error analysis that separates sampling, learning and time-stepping errors: up to the learning error, NEOPIC has the error structure of standard Monte Carlo particle methods. We show that the particle noise in the training data acts as an implicit regularizer, which gives guidance on the number of training particles in relation to the particle numbers used at inference. A complexity analysis shows that NEOPIC is an $O(N_p)$ method whose constants are set by the architecture alone, independently of any mesh or tree. Numerical experiments on electrostatic benchmarks in two and three dimensions, including Landau damping, two-stream and bump-on-tail instabilities, a cyclotron mini-app and a Penning trap, confirm the theoretical findings. A single model trained on several mini-apps generated by different reference methods captures all of them and remains stable up to twice the training time window, which points towards foundation models for fields in particle-based methods.
This article presents two new multilevel hierarchical matrices to perform fast kernel matrix-vector products in $d$ dimensions $(d>1)$, namely $\mathcal{H}_{hw}^{\prime}$, $\mathcal{H}_{hw}^{\prime \prime}$. We introduce a novel hybrid weak admissibility condition to construct these hierarchical matrices. More precisely, at coarser levels of the tree (i.e., larger blocks), we use a weak admissibility condition, wherein \textbf{far-field} and \textbf{vertex-sharing} clusters are considered admissible, while at the finer levels (i.e., smaller blocks), the \textbf{adjacent} clusters are treated as admissible. The admissible blocks at coarser levels of $\mathcal{H}_{hw}^{\prime}$ and $\mathcal{H}_{hw}^{\prime \prime}$ are represented using semi-nested and fully nested bases, respectively. In both hierarchical matrices, the finer levels at which adjacent clusters are treated as admissible are constructed in a non-nested fashion. Thus, they can be viewed as a hybrid of $\mathcal{H}$ and $\mathcal{H}^2$ matrices. We use adaptive cross approximation (ACA) and nested cross approximation (NCA) to construct non-nested and nested bases, respectively. The fast matrix-vector products are performed using the proposed hierarchical matrix formats. To evaluate the performance and usefulness of the proposed hierarchical matrices, we conduct some numerical experiments. The numerical results show that they outperform their uniform admissibility counterparts in terms of storage and matrix-vector product cost, demonstrating the effectiveness of the proposed hybrid formats. We also show that the proposed formats are competitive with the standard (or strong) admissibility-based $\mathcal{H}^2$-matrices. Furthermore, we showcase some interesting applications of the proposed hierarchical matrices.
We study multigrade deep learning (MGDL) as a principled framework for structured error refinement in deep neural networks. While the approximation power of neural networks is now relatively well understood, training very deep architectures remains challenging due to highly nonconvex and often ill-conditioned optimization landscapes. In contrast, for relatively shallow networks, most notably certain one-hidden-layer ReLU models, training admits convex reformulations with global guarantees under appropriate settings, motivating learning paradigms that improve stability while scaling to depth. MGDL builds on this insight by training deep networks grade by grade: previously learned grades are frozen, and each newly added grade-wise subnetwork is composed on top of the previously learned grades and trained to fit the residual left by the current approximation, yielding a structured and interpretable hierarchical refinement process. We develop an operator-theoretic foundation for MGDL and prove that, for any continuous target function defined on a hypercube, there exists a fixed-width multigrade ReLU scheme whose residuals are pointwise nonincreasing in magnitude and converge uniformly to zero, with strict $L^p$-norm decay at every nontrivial grade for $p\in [1,\infty)$. To the best of our knowledge, this work provides the first rigorous constructive approximation guarantee showing that a grade-wise residual refinement scheme can achieve vanishing error in a fixed-width multigrade ReLU architecture.
We present an adaptive bitwidth allocation method for quantizing low-rank factors in hierarchical matrices in order to decrease storage while meeting a given constraint on their reconstruction error. For the bitwidth allocation we formulate a discrete storage minimization problem and show that an optimal assignment can give larger singular values at least as many bits as smaller ones. We further develop a greedy algorithm that meets the specified max-norm quantization error budget. We investigate a hierarchical block partition in which all of the non-diagonal blocks are represented by quantized low-rank factors and only dense blocks remain dense; in this way a global reconstruction bound related to $η_k$ and $R_k$ is obtained, with the former and latter being bounds on the truncation and reconstruction arithmetic errors, respectively. Extensive simulations, carried out with kernel and BEM matrices under block-relative budgets, indicate that the proposed method achieves a reduction in storage of up to $44.3\times$ when compared to storing the complete matrix densely in FP64, and a reduction of up to $5.70\times$ when compared to the FP64 low-rank factors within the same partition, together with measured reconstruction errors that remained below the maximum norm bounds in all the hierarchical tests.
Ana Čukarska, Maximilian Gollwitzer, Viola Muning Li, Vladyslav Fediukov, Zahra Monfared, Felix Dietrich
Learning continuous-time dynamical systems from data is an active area of research. The Koopman formalism is often used here, because it is a flexible approach that is informative about the system's fundamental properties. The formalism represents dynamical systems through linear operators acting on a suitable function space. A practical difficulty arising in the numerical approximation of Koopman operators is the choice of a ``dictionary'', meaning an effective finite-dimensional representation of the function space. In this paper, we focus on learning dynamical systems in continuous-time, and so we consider numerical approximations of the Koopman generator instead of the Koopman operator. Our main contribution is a data-driven random-feature method to obtain a parsimonious and effective function representation for Galerkin approximation methods such as gEDMD. The representation we construct is data-informed, which allows us to improve upon data-agnostic methods. We empirically demonstrate the suitability of our method for forecasting tasks and system analysis. The numerical experiments show that our approach is accurate, and faster than iterative optimization. Further, we prove convergence of our approximation of the vector field and the resulting trajectories. In addition to the numerical results, we also make theoretical contributions, including a general convergence criterion for the semigroup generated by the gEDMD approximation and a counterexample showing that gEDMD does not generally preserve exponential boundedness. With our contribution, we advance the frontier of efficiently learning continuous-time dynamics from data, towards interpretable and mathematically robust modelling based on the Koopman framework.
Jan Blechschmidt, Oliver G. Ernst, Moritz Poguntke, Björn Sprungk
Sampling from the posterior is the central task of computational Bayesian inverse problems. The standard workhorse in Bayesian inference---Markov chain Monte Carlo (MCMC)---is sequential, yields correlated samples, and must be rerun for each new observation. Conditional flow matching offers an amortized alternative: a transport map, trained once on joint samples of parameter and data, that yields independent approximate posterior samples for any observation at negligible online cost, without new likelihood evaluations. We give a careful, MCMC-literate assessment of flow matching for PDE-based inverse problems with function-valued parameters. Exploiting the flow's tractable density, we derive computable accuracy estimates for the learned posterior in total variation distance and Kullback--Leibler divergence and, moreover, propose a hybrid sampler that is made asymptotically exact by Metropolization. We validate the accuracy estimates and demonstrate amortization in several numerical examples, including electrical impedance tomography and a likelihood-free Lotka--Volterra model. Flow matching reproduces posterior means and credible bands closely at a fraction of the forward solves required by MCMC, but its accuracy degrades for concentrated posteriors and for observations that are atypical under the prior.
Reduced-precision arithmetic offers significant opportunities to improve performance, memory usage, and energy efficiency in numerical applications, provided that numerical accuracy is preserved. This work investigates automated precision tuning through customized floating-point formats with user-defined exponent and significand sizes, enabling the emulation of emerging low-precision formats and the exploration of non-standard precision configurations within a unified mixed-precision framework. The proposed methodology, implemented in the PROMISE precision autotuning tool, combines numerical validation with a systematic search to generate program variants that satisfy user-defined accuracy requirements. To address the computational cost of this exploration, a containerized benchmarking framework supports parallel execution across multiple algorithms and parameter configurations. The approach is evaluated on a suite of numerical programs, including linear solvers and applications from the Rodinia benchmark. Results show that a substantial proportion of variables can be safely reduced to lower precision while preserving accuracy, indicating that standard double precision is often over-provisioned. These findings highlight the potential of automated precision tuning to derive efficient mixed-precision configurations tailored to application-specific accuracy requirements.
Low-rank methods can reduce the memory and computational requirements of deep neural network training in approaches such as GaLore. Randomized range finders offer an attractive alternative to singular value decompositions, particularly when the required rank is determined adaptively from a prescribed approximation tolerance. We introduce a blocked adaptive randomized range finder based on a Householder QR decomposition and an optimized hybrid CPU--GPU implementation. Numerical experiments show that it preserves orthogonality and approximation accuracy in regimes where block Gram--Schmidt without reorthogonalization becomes unstable. The blocked formulation exposes matrix--matrix operations and enables overlap of CPU panel factorization with GPU updates. On an NVIDIA GH200, the overlapped implementation reduces the runtime for the largest tested matrix from 9.91 s on the CPU to 0.407 s. The method provides a stable and efficient building block for low-rank approximation on heterogeneous systems with applications in computational science and engineering.
Instantaneous Quantum Polynomial-time (IQP) circuits are promising candidates for near-term quantum computing due to the conjectured classical hardness of their sampling task and their shallow, highly parallelizable structure. Despite extensive studies of their complexity-theoretic properties, the role of circuit connectivity in their trainability as variational ansatze remains poorly understood. In this work, we present an investigation of IQP circuits for Hamiltonian optimization that traces a mechanistic chain linking circuit connectivity, entanglement generation, and gradient scaling: increasing connectivity drives the circuit output toward Haar-like statistics and volume-law entanglement, and this increase in entanglement is consistently associated with an exponential decay of the gradient variance and degraded optimization performance, while sparser architectures retain stable gradients at the cost of global expressibility and entanglement capacity. This connectivity-entanglement-gradient mechanism produces a nontrivial connectivity-trainability trade-off, demonstrating that effective IQP-based Hamiltonian optimization requires a balance between global expressibility and trainability rather than maximal expressivity alone. In particular, on our benchmarks, a circuit generating no entanglement already matches the optimization performance of more highly connected circuits, and a circular topology retains entanglement capacity while preserving trainability, making it a broadly effective choice among the architectures considered.
Direct evaluation of the defining product of the Cayley-Dickson algebra $A_n$ of rank $N:=2^n$ over a unital commutative ring $R$ in which $2$ is invertible requires a quadratic number of arithmetic operations. This paper gives a uniform bilinear algorithm for multiplication in this algebra costing $O(N \log N)$ arithmetic operations in $R$ and $O(N)$ auxiliary storage when the doubling parameters of $A_n$ are units whose square roots lie in a fixed commutative coefficient algebra $S$ finitely generated as an $R$-module and admitting an $R$-linear map $π\colon S \to R$ with $π(1) = 1$. In particular, over $\mathbb{R}$ and every finite field $\mathbb{F}_{q}$ of odd characteristic, a quadratic field extension contains a square root of every possible doubling parameter, so the algorithm applies to every Cayley-Dickson tower over these fields, and the tensor rank of multiplication in $A_n$ over them is at most $(9n-21)2^n+28$ when $n\ge2$. When the top doubling parameter is a non-square, as in the standard real tower, the bound is $(9n-15)2^{n-1}+10$. To the author's knowledge, the algorithm presented in this paper is the first uniform subquadratic algorithm for multiplication in $A_n$. When $R=\mathbb{R}$ with every doubling parameter equal to $-1$, the real quasilinear schedule uses fewer total arithmetic operations than both direct multiplication and a uniform algorithm of Cariow-Cariowa for $n\ge4$. An implementation for the standard real tower is available in the MIT-licensed C11 library fastCD, with a NumPy-backed Python interface. In single-core benchmarks against direct multiplication and the uniform Cariow-Cariowa method, the quasilinear implementation had the lowest mean time of the three at every tested dimension $N\ge32$, for both single and batched products, and was over $17$ times faster than direct multiplication at $N=1024$.
Enforcing prescribed global integral constraints in mesh-free neural PDE solvers is challenging in high-dimensional domains. Existing projection methods for spatial integrals are often tied to fixed grids or uniform quadrature, which can conflict with randomly sampled physics-informed neural networks (PINNs) and scale poorly with dimension. High-order differential operators also increase reverse-mode automatic differentiation memory costs. We propose Stochastic Dimension Implicit Functional Projection (SDIFP), a quadrature-level framework for enforcing prescribed first and second spatial moments. SDIFP replaces tensor-product nodal projection by a global affine correction of the neural-network output, with two scalar coefficients determined from a weighted quadrature rule. Under positive target variance and nonzero empirical raw variance, this correction is the nearest-point projection, in the weighted quadrature norm, onto the empirical two-moment constraint set. Thus, the prescribed moments are exact for the selected quadrature rule, while continuum errors are quadrature errors of the corrected field. For decomposable high-dimensional linear operators, SDIFP combines affine moment correction with stochastic operator-subset sampling. With independent residual and derivative sampling and conditionally unbiased coefficient-gradient estimation, the resulting estimator is unbiased for the specified quadrature-based residual objective; the shared-subset fast mode is biased in general. SDIFP avoids tensor-product quadrature for moment enforcement, separates forward quadrature evaluation from the reverse-mode graph, and retains pointwise inference efficiency once the affine coefficients are fixed or precomputed.
Oscar Mickelin, Yunpeng Shi, Amit Singer, Eric J. Verbeke
We develop a fast method for the rigid alignment of voxelized three-dimensional density functions over $\mathsf{SO}(3)$ using sliced optimal transport. Our approach efficiently transforms each density map into a tensor of one-dimensional quantile functions, taken from central slices over the unit sphere, such that the Euclidean distance of the transformed volumes approximates the sliced 2-Wasserstein distance. Our formulation is therefore compatible with existing, fast, voxel-based Euclidean alignment algorithms. We provide theoretical analysis and numerical experiments demonstrating that the sliced 2-Wasserstein distance is a robust metric for comparing three-dimensional shapes and for computing their rigid alignment.
In variable-step operator learning, agreement at observed step sizes does not by itself determine predictions at other temporal resolutions. We propose an operator whose step-size dependence is prescribed by a first-order Duhamel approximation of semilinear dynamics. It combines exact linear propagation with a local nonlinear function of the field and its derivatives, allowing evaluation below the observation interval without supervision densely sampled across step sizes. Under suitable conditions and at fixed spatial resolution, we prove that agreement between models at two step sizes controls their prediction differences throughout a prescribed range. We also give conditions for local identifiability of both components from two step sizes, supporting interpretation of the learned structure. A latent extension applies the same construction to learned fields. We evaluate the base model on two synthetic problems and the extension on two real datasets against several baselines.
We present new representations of Gauss--Legendre polynomials and their derivatives in the shifted power basis and in bases related to symmetric orthogonal Jacobi polynomials. Using these representations and certain recurrence relations, we propose $O(n^2+dn)$ methods for evaluating a Gauss--Legendre curve of degree $n$ in $\mathbb E^d$. We also propose algorithms for multipoint evaluation with computational complexity $O(Mdn+dn^2)$, where $M$ is the number of evaluation points. Numerical tests show that the approach using the shifted power basis is very fast, but it works well only for small $n$ (approximately up to 20), whereas the method employing Jacobi polynomials gives good numerical results even for large $n$ (around 100).
It is well-known that a $d$-dimensional polynomial Bézier curve of degree $n$ can be subdivided into two segments using the famous de Casteljau algorithm in $O(dn^2)$ time. We show that it is possible to do this in $O(dn\log{n})$ time using the fast Fourier transform and its inverse. The direct application of the method is numerically unstable for larger values of $n$. Therefore, we also introduce a modified version based on parameter scaling, which retains the same asymptotic computational complexity. This scaling approach can improve numerical accuracy in some cases (for $n\lesssim 50$), as demonstrated by our experiments. Moreover, the new method has a nice property: if a Bézier curve is extended by an additional control point, the subdivision can be updated in $O(d)$ time. The method can be adapted to subdivide rational Bézier curves and rectangular Bézier surfaces, as well as to compute the derivatives of Bézier curves. Although motivated by applications in CAGD, the proposed methods may also be useful in other areas involving polynomial reparameterization.
Physics-informed holomorphic neural networks (PIHNNs) exploit complex-analytic representations to satisfy homogeneous PDEs by construction, but their extension to nonhomogeneous problems is nontrivial. We propose a hybrid physics-informed holomorphic neural network (HPIHNN) framework based on a decomposition of the solution into particular and homogeneous components. The particular component is computed by a spectral method on a simple extended domain, while the homogeneous component is approximated by a PIHNN trained only on the boundary with modified boundary data. The resulting formulation avoids meshing complex geometries and retains the dimensionality-reduction advantage of holomorphic representations. The method is applied to nonhomogeneous Poisson, Stokes, biharmonic, linear elasticity, and Helmholtz equations. Numerical experiments covering complex geometries, corner singularities, and nearly incompressible regimes demonstrate that the proposed method accurately resolves the corresponding solutions and, in the considered tests, generally achieves lower errors and shorter training times than standard PINNs. The results indicate that coupling spectral particular solutions with holomorphic neural approximations provides an efficient approach to nonhomogeneous PDEs on complex domains.
Kernel methods face two recurrent bottlenecks: the cubic cost of inverting an $n\times n$ kernel matrix, as in kernel ridge regression, and the evaluation of kernel integrals, as required in kernel quadrature. Existing low-rank kernel and Gaussian process approximations primarily address the former bottleneck, but do not control the downstream quadrature error when both the kernel matrix and kernel mean are approximated, particularly when the integrand lies outside the kernel reproducing kernel Hilbert space. We address both challenges through an operator-adapted finite-feature approximation. Specifically, we construct a finite-dimensional working kernel from the eigenvalues and eigenfunctions of a second-order differential operator. We establish convergence guarantees for the resulting finite-feature estimators in kernel quadrature, with extension to other related kernel methods. Our results show that, with $m \lesssim n$ basis functions, the truncated estimator attains the corresponding benchmark statistical rates. The theory covers fixed/random designs, noisy/noiseless observations, oversmoothing/undersmoothing regimes, and Matérn/Gaussian-type weights. We provide explicit error bounds that quantify the roles of the spectral cutoff, regularization parameter, lengthscale, sampling geometry, and operator boundary condition. We verify the convergence results in numerical experiments.
Xinyuan Li, Matthew Zahr, Panos Stinis, Zecheng Zhang
Multi-fidelity learning uses low-fidelity information (coarse-scale PDE models) together with limited high-fidelity data (fine-scale solutions) to obtain high-accuracy PDE solutions. It is a data-driven approach to downscale reduced-order models to solve PDEs efficiently and accurately. Standard multi-fidelity learning focuses on one operator or PDE; given a new PDE, retraining using additional coarse-fine solution pairs may be required, which is costly. In this work, we propose an in-context multi-fidelity learning framework for learning multiple multi-fidelity operators with one shared model. The model is trained using paired low- and high-fidelity examples as demonstrations from multiple operators, together with a new coarse-scale solution query, to predict the corresponding fine-scale solution. Once trained, given a new PDE operator, a few (fewer than five) coarse-fine demonstration pairs, and a low-cost coarse query, the model predicts the corresponding fine-scale solution without fine-tuning. When the target operator appears in the training dataset but has a different prediction goal, such as predicting beyond the training time interval, no additional demonstrations are needed. Instead, the model uses training coarse-fine pairs from that PDE and information from other training operators. Testing then only requires generating a low-cost coarse-scale solution query, without fine-tuning. Notably, our smallest model contains only 255,364 trainable parameters, corresponding to approximately 456 parameters per training operator for 560 operators.
David Voderholzer, Leonard Vetter, Thomas O'Leary-Roseberry
Parametric PDEs in science and engineering are often posed on irregular meshes with geometries that may contain holes or complex boundaries, yet popular operator-surrogate architectures such as the Fourier Neural Operator (FNO) rely on structured geometries for fast computational kernels. We introduce Adaptive Modal Operator Refinement (AMOR), which replaces FNO's Fourier-based spectral convolutions with reduced-basis blocks and supports both data-driven and prescribed input/output bases through a common encode/decode interface. AMOR thus applies unmodified to arbitrary meshes and, because its spectral parameterization has no tensor-product dependence on the spatial dimension, it scales favorably to three dimensions, where FNO's spectral weights grow as $m^3$ with the number of spectral modes $m$. We prove that AMOR is a universal approximator of continuous operators. On two-dimensional structured-grid benchmarks, AMOR is competitive with FNO overall and surpasses PCANet: it outperforms FNO on nonlinear diffusion and Helmholtz problems, while a static snapshot basis keeps it from matching FNO's accuracy on the time-dependent Kolmogorov flow problem. On three-dimensional problems it matches or exceeds FNO's accuracy at lower inference cost on linear elastostatics while trailing FNO on a turbulent lid-driven-cavity flow and scales favorably to problems demanding high capacity. On irregular meshes, where FNO is not directly applicable, AMOR reduces the relative $L^2$ error over the FNO extension Geo-FNO by $31\%$ (airfoil) and $42\%$ (elasticity). Finally, AMOR maintains strong performance across three physically distinct, experimentally motivated loading regimes on a topologically nontrivial three-dimensional aircraft-rib mesh, demonstrating its applicability to realistic engineering problems without modification.
Smoothed particle hydrodynamics (SPH) is well suited to large-deformation fluid problems, and the growing use of automatic differentiation in engineering applications has motivated interest in differentiable SPH solvers. In this study, a composable differentiable SPH framework, CraftSPH, has been developed. In CraftSPH, major computational operations are implemented as independent modules that users can combine to construct solvers for different purposes. This design supports high-accuracy discretization and both explicit and implicit computation within a common framework, while allowing individual components to be maintained and extended. Selected numerical operations also support automatic differentiation (AD) for physical-parameter estimation and potential coupling with deep learning models. These capabilities were demonstrated through forward and inverse analyses. Forward analyses using an implicit SPH scheme with high-accuracy discretization showed good agreement with analytical and reference solutions for Poiseuille flow and a rising bubble, while capturing the key experimental trends in two- and three-dimensional dam-break problems. In inverse analyses, explicit and implicit SPH schemes were used to estimate kinematic viscosity in Taylor-Green vortex flow and Reynolds number in lid-driven cavity flow. The results also showed that the choice of operations included in the computational graph affects the computed gradients. Together, these findings demonstrate that CraftSPH enables users to assemble problem-specific solvers for forward and inverse analyses within a common framework.
Alba Muixí, Lluís Monforte, Antonio Rodríguez-Ferran, Sergio Zlotnik
This work proposes a hyper-reduced order model for coupled hydro-mechanical simulations of elastoplastic deformation and water flow in partially saturated soils. These simulations are used to monitor and assess the safety of geotechnical structures such as earthfill dams, where calibration, uncertainty quantification and scenario analysis require solving the same transient model many times for different parameter values. The displacement and the pore pressure are approximated in a field-split reduced basis obtained by proper orthogonal decomposition, and the Galerkin-projected equations are integrated with decomposed cubature rules, one per balance equation, built with the energy-conserving sampling and weighting (ECSW) method. Each residual is split into an affine part and a nonlinear correction: the affine part is computed exactly, at negligible online cost, so that the cubature rule only needs to approximate the correction, which is localized in space. The construction is independent of the constitutive model choices and is implemented in FEniCSx by assembling the original variational forms on the small set of elements selected by ECSW. The method is tested on two- and three-dimensional dams retaining tailings and water, simulating a raise of the dam and the subsequent consolidation phase, with the friction angle, stiffness and hydraulic conductivity of the dam material and the height of the raise as varying parameters. The hyper-reduced model reproduces the displacement and pressure fields with errors below $0.5\%$ with respect to the finite element solution, achieving speed-ups of up to $43\times$ in two dimensions and $77\times$ in three. These performance gains grow with the extent of the plastic yielding, and the selected elements concentrate in the plastified and unsaturated zones, whose size is set by the physics rather than by the mesh resolution.
Walk-on-Spheres (WoS) enables mesh-free solutions of elliptic PDEs, but conventional methods repeatedly discard sampled trajectories, limiting their efficiency in interactive applications. We introduce reusable transport skeletons, a persistent representation that organizes stochastic walks into a shared computational graph. Our framework enables efficient reuse across spatial queries and boundary conditions, while supporting progressive refinement, structure-preserving compression, and local adaptation to dynamic geometry. It further enables interactive editing of transport between selected regions. Experiments demonstrate substantial computational savings while preserving solution accuracy. By transforming transient stochastic trajectories into reusable transport structures, our method provides a flexible foundation for interactive Monte Carlo PDE solving.
Changjian Zhou, Junfeng Fang, Negin Yousefpour, Peng Wu, Bin Yan, Guillermo A. Narsilio
Neural operator models trained on simulation data often lose accuracy when applied to experimental measurements due to the sim-to-real gap. Standard fine-tuning with limited real data can reduce this gap, but it may also damage the core physics-relevant representations learned during pretraining. Although knowledge-preserving adaptation has been widely investigated in vision or language tasks, it remains unclear whether these methods are suitable for neural operators whose architectures and protected knowledge are fundamentally different. Neural operators need to preserve core-scale physical structures rather than semantic or visual features. We propose PhysGuard, a physics-preserving framework for accurate sim-to-real adaptation of neural operators. Specifically, PhysGuard uses the empirical Fisher Information Matrix computed on simulation data to identify physics-critical parameter directions, then restricts fine-tuning updates to directions that do not interfere with them. A layer-wise Gram-matrix formulation makes this efficient for models with millions of parameters, while an adaptive threshold automatically determines the protected subspace size. A spectral probe experiment shows that the dominant Fisher directions are strongly associated with low-frequency output structures. Experiments on benchmark across four neural operator architectures and different physical systems show that PhysGuard performs strongly on most evaluation metrics compared to baselines. The benefits are most evident under severe domain shift, where it reduces low-frequency error by up to 32% compared to standard fine-tuning while maintaining adaptability. Our code is available at https://github.com/Data-Driven-Computational-Geotechnics/PhysGuard.
Stationary tail computation for high-dimensional reflected Brownian motions (RBMs) remains computationally challenging. Building on NN4MGF's shared-coordinate basic-adjoint-relationship (BAR) framework, we develop Holomorphic Cauchy Networks (HCNs). Native complex-valued rational Cauchy subnetworks enforce holomorphicity on pole-free domains, removing sampled Cauchy-Riemann penalties and their input-derivative graphs. Exact origin anchoring and boundary-coordinate exclusion enforce structural identities, while anchored interactions represent dependent transforms. Under holomorphic-log assumptions, we establish simultaneous value and first-derivative approximation for the constrained architectures and a mixed-derivative lower bound for coordinate-only representations under nonzero cross-covariance. Closed-form real-axis derivatives and sequential backpropagation further control memory. On the same NVIDIA H100 NVL, the tested 30D HCN configuration uses approximately tenfold less peak allocated memory and has 9.7-fold lower mean update time than NN4MGF; extrapolation to the respective scheduled update counts gives a projected full-schedule update-time ratio of 58.2. With 64 million sampled base inputs per high-dimensional run, HCN reaches 50 dimensions on a laptop GPU with 6.367% mean relative error over thresholds 1 through 8. Degree-matched comparisons yield lower reported mean relative errors over the primary evaluation ranges, including 1.422% versus 6.684% on the 2D Harrison non-product benchmark. An exploratory 10D correlated example illustrates the interaction extension.
Valentin Leplat, Daniil Merkulov, Aleksandr Katrutsa, Daniel Bershatsky, Olga Tsymboi, Ivan Oseledets
Stochastic first-order methods are widely used for large-scale optimization and machine learning, where their performance often depends on the learning rate. We study NAG-GS, a stochastic momentum method based on a semi-implicit Gauss--Seidel discretization of a Nesterov-type flow. For strongly convex quadratics, we determine the exact range of stable step sizes for fixed damping and the step size that minimizes the worst-case spectral radius. When the damping equals the strong-convexity parameter, the method recovers the accelerated asymptotic factor $(\sqrtκ-1)/(\sqrtκ+1)$, where $κ$ is the condition number. We also prove that the quadratic iteration converges linearly with dynamic damping when its limiting fixed-damping iteration is stable. For general strongly convex, possibly nonsmooth objectives, we establish an $O(1/k)$ bound on the expected objective gap of the last iterate, using a conditionally unbiased stochastic first-order oracle with bounded second moment. The result holds with dynamic damping from any $γ_0>0$. We show that NAG-GS can also be written as a Heavy-Ball iteration with explicitly determined step sizes and momentum coefficients. We test the stochastic result on nonsmooth hinge-loss classification problems, and further evaluate NAG-GS on quadratic problems, logistic regression, convolutional neural networks, RoBERTa, and Vision Transformers. The experiments show that NAG-GS can remain effective at learning rates where standard optimizers become unstable or lose accuracy.
Nearly incompressible materials arise in soft tissues, elastomers, gels, and other soft solids. Yet, numerical modeling remains challenging because displacement-based formulations can exhibit volumetric locking near incompressibility. Mixed displacement-pressure formulations address this difficulty but require compatible approximation spaces satisfying the inf-sup condition. In mixed formulations, discontinuous pressure spaces are attractive since pressure degrees of freedom can be condensed locally. In high-order isogeometric analysis, however, their Galerkin stability requires pressure-grid coarsening that increases with polynomial degree. Here, we develop a variational multiscale (VMS) stabilized mixed isogeometric formulation that retains this locality while substantially relaxing the restriction. The key observation is that VMS stabilization controls nonconstant elementwise pressure modes and is unconditionally stable with respect to its parameter, provided the elementwise constant pressure component satisfies the inf-sup condition. Consequently, for maximally continuous displacement splines, a pressure grid coarsened by a factor of two satisfies the constant-mode inf-sup criterion independently of polynomial degree. Local condensation then yields a displacement-only system, facilitating use in existing displacement solvers. The formulation is further extended to immersed isogeometric analysis with extended B-splines, pressure-cell aggregation, subdivision-based cut-cell quadrature, and weak enforcement of essential boundary conditions. Discrete inf-sup studies and convergence analyses in the Kirsch and Cook's membrane benchmarks demonstrate stability and show that the immersed treatment preserves mixed stability properties. The proposed framework therefore provides a local, high-order, and robust approach for nearly incompressible problems in mesh-conforming and immersed settings.
The curse of dimensionality associated with solving the Hamilton-Jacobi-Bellman equation for the value function remains a key obstacle for real-time optimal feedback control of high-dimensional nonlinear dynamical systems. This is especially severe in aerospace applications such as relaminarizing unsteady turbulent flows and recovering from aerodynamic stall where nonlinear dynamics play a significant role. We introduce Balanced Active Subspaces for Optimal Control (BASOC), a model reduction framework which seeks to find directions in state space along which the value function is most sensitive, and the distribution of states exhibits maximum variance along optimal trajectories. BASOC leverages trajectory optimization and the Pontryagin minimum principle to obtain state and co-state (value function gradient) data along locally optimal trajectories for constructing state and gradient covariance matrices. Balancing these matrices provides an oblique projection for approximating the value function, which can also be used to construct reduced-order models and surrogate feedback control laws. We show that for linear-quadratic optimal control problems, BASOC becomes equivalent to a snapshot method for initial condition balancing of the optimally controlled closed-loop realization. The framework is demonstrated on a nonlinear three-dimensional system where POD and balanced truncation are known to perform poorly, and on the nonlinear complex Ginzburg-Landau equation to demonstrate feasibility for larger-scale systems.
We propose an optimization-informed deep neural network approach, named iUzawa-Net, aiming for the first solver that enables real-time solutions for a class of nonsmooth optimal control problems of linear partial differential equations (PDEs). The iUzawa-Net unrolls an inexact Uzawa method for saddle point problems, replacing classical preconditioners and PDE solvers with specifically designed learnable neural networks. We prove universal approximation properties and establish the asymptotic $\varepsilon$-optimality for the iUzawa-Net, and validate its promising numerical efficiency through nonsmooth elliptic and parabolic optimal control problems. Our techniques offer a versatile framework for designing and analyzing various optimization-informed deep learning approaches to optimal control and other PDE-constrained optimization problems. The proposed learning-to-control approach synergizes model-based optimization algorithms and data-driven deep learning techniques, inheriting the merits of both methodologies.
Anastasis Kratsios, Tin Sum Cheng, Bum Jun Kim, Daniel Roy
Transformers often transfer across related tasks with little fine-tuning; we study this question in noisy nonparametric regression for continuous maps $f,g:[0,1]^d\to\mathbb{R}^D$, representing the source and shifted tasks, using a single-query, single-head transformer $\mathcal{T}(\cdot)=V\operatorname{Softmax}(-α(\cdot,K))$ with $\ell^{\infty}$-based alignment scores $α$. Rather than assume that gradient-based training reaches an empirical risk minimizer, we give a model-specific two-stage procedure \textbf{certifying} source-task recovery and out-of-distribution (OOD) adaptation. With high probability, the resulting attention head recovers $f$ to arbitrary accuracy $\varepsilon>0$ in every fixed $L^p([0,1]^d)$ norm, $1\le p<\infty$, from noisy observations. For Lipschitz targets, its rate is minimax-optimal up to logarithmic factors without assuming oracle access to an asymptotic ERM. If $D=1$ and $f$ is Pfaffian than the values matrix factors into $AC$ with $A$ scaling at-most linearly $\mathcal{O}(\varepsilon^{-1})$ in the reciprocal approximation error $\varepsilon>0$. Importantly, our training procedure encodes all our identified (potentially) shared geometric structure from $f$ only into the values matrix factor $C$ and the keys matrix $K$ allowing for rapid access, and all $f$-specific structure only into $A$. Our main result concerns the OOD phase; where we identify geometric structures, which if shared by the source $f$ and OOD $g$ task; then only only the values matrix factor $A$ needs to be fine-tuned to guarantee an $L^p$-recovery of $g$. Critically, fine-tuning only $A$ is a convex problem which involving an \textit{at-most linear} number of parameters $\mathcal{O}(\varepsilon^{-1})$; dramatically cutting down on the worst-case scaling of $\mathcal{O}(\varepsilon^{-\mathbf{d}})$ unavoidable in the source-training phase.
Abhinab Bhattacharjee, Andrey A. Popov, Amit N. Subrahmanya, Adrian Sandu
This work proposes a new Multi-Fidelity Particle Flow Filter (MFPFF) based on a linear control variate framework and multilevel Monte Carlo methods. Conventional particle filters often suffer from weight collapse and particle degeneracy in high-dimensional spaces. Particle flow filters address this by smoothly transporting prior particles toward the posterior distribution, thereby avoiding the need for resampling. This work constructs multi-fidelity particle flow filters by augmenting a small number of expensive high-fidelity particles (full-order model runs) with a large number of inexpensive low-resolution particles (reduced-order surrogates constructed via intrusive Galerkin projection). These surrogates act as control variates to enhance the statistical estimate of the full-order posterior, yielding enhanced analysis accuracy at comparable computational cost. The work focuses on two particle flows based on log homotopy, the stochastic Gromov flow and deterministic Exact flow. Numerical experiments with the Korteweg-De Vries (KdV) and Quasi-Geostrophic (QG) test problems demonstrate that the proposed MFPFF achieves accuracy comparable to high-fidelity simulations while significantly reducing computational expense.