Papers

This paper addresses infinite-dimensional Bayesian inference for inverse problem of partial differential equations with model parameters in infinite-dimensional Hilbert space. To effectively incorporate prior information, we propose a novel continuous normalizing flows based infinite-dimensional model. Specifically, by introducing a well-defined neural ordinary differential equation in infinite-dimensional space, a simple reference measure can be transformed into a more complex measure which encodes the prior information. A corresponding theoretical framework is established to ensure the well-posedness of our proposed Bayesian prior in infinite-dimensional space. We also provide training methods of the prior for two distinct data settings, along with two sampling algorithms for the resulting Bayesian posterior. The proposed framework is applied to three representative inverse problems: the simple smooth inverse problem, inverse scattering problem, and the inverse heat conduction problem. Numerical experiments support the theoretical analysis and demonstrate the efficiency of the proposed algorithms.

Affine-invariant ensemble samplers are widely used in Bayesian applications. However, their quantitative convergence theory, in particular geometric ergodicity, remains a basic open question. We study the affine invariant ensemble Langevin dynamics, an interacting particle system that uses the empirical covariance of the whole ensemble as a preconditioner. While effective in practice, theoretical understanding of this method is not available beyond plain qualitative convergence in total variation, with a central difficulty being that the empirical covariance can become degenerate. For potentials with bounded Hessian that are strongly convex outside a ball, we prove geometric ergodicity using a novel Lyapunov function that combines the inverse covariance with a coercive exponential energy. We then show that directly applying the Euler–Maruyama scheme can diverge with positive probability, even for a one-dimensional Gaussian target. This motivates a covariance-trace regularization in continuous time. We prove geometric ergodicity of the regularized diffusion and, for sufficiently small step size, of its unadjusted Euler–Maruyama discretization. We also show that the invariant distributions of the discretization converge weakly to the product target distribution as the step size tends to zero. Lastly, we provide an explicit discretization that is both affine-invariant and geometrically ergodic, based on an adaptive step size clamping.

This paper introduces a novel, efficient class of parallel direct and indirect iterative schemes to solve general obstacle and free boundary value problems. The uniqueness of the solution for the direct parallel method is established under the assumption that the model problem yields an $M$-matrix. The convergence analysis of the indirect approach is predicated on minimizing the discrete energy functional at each iterative stage of the linear approximation process. Under the aforementioned setting, we establish theoretical convergence results for both smooth and nonsmooth energy functionals defined over a convex set, in the sense of Ferris and Mangasarian \cite{ferris1994parallel}. For the numerical computation of the one-dimensional obstacle problem, we adopt a direct generalized parallel approach based on the SPIKE algorithm \cite{polizzi2006parallel}. To accelerate convergence and reduce computational complexity, a fast recursive version of the scheme is generalized to the nonsmooth case. For two- and three-dimensional obstacle problems, we employ a directional splitting method that treats each directional subproblem as a one-dimensional obstacle minimization problem. By framing the energy functional minimization as a parabolic time-dependent problem and utilizing an Armijo time-stepping rule during the solution update process, we successfully obtain results for higher-dimensional obstacles. Additionally, this study explores the possibility of extending the algorithm into a constrained quadratic programming optimization solver. We also evaluate the framework on image deblurring phenomena under the aforementioned setting, successfully recovering the original images. Finally, numerical illustrations are provided to validate the theoretical results.

Many biomedical challenges can be posed as tensor completion problems where the observed entries of a multidimensional array (a tensor) are used to impute the missing values. In such settings, incorporating side information about the modes of the tensor, such as gene-gene similarity, can significantly enhance the solutions of the completion problem. Most existing tensor completion methods can only incorporate side information in the form of matrices. In this study, we introduce a novel framework to incorporate side information in the form of tensors. Our new approach, called Coupled Tensor-Tensor Completion (CTTC), leverages the hidden connections among multimodal tensors to improve tensor completion performance. In addition to practical utility, CTTC has theoretical foundations in distance metric learning and group theory. We derive an alternating algorithm to solve the CTTC optimization problem and establish its convergence to a stationary point. Finally, we show that CTTC outperforms state-of-the-art tensor completion methods at predicting drug effects. Results: Compared with other tensor completion methods, including HaLRTC, CTRC, Cell, and NTDDR, CTTC demonstrates superior run-time and RSE tensor completion accuracy on two benchmark datasets, DTD and LINCS.

We present a high-order multiphysics solver for the simulation of microwave-heated flows. The solver couples a discontinuous Petrov-Galerkin (DPG) finite element method for the time-harmonic Maxwell equations with continuous Galerkin finite element methods for the heat equation and the incompressible Navier-Stokes equations. We validate the electromagnetic solver against multiple benchmark problems: wave propagation in a rectangular waveguide, a cavity problem with a singular solution, and a microwave-heated obstacle problem, comparing our results against numerical and experimental data from the literature. The results confirm the validity of the implementation and demonstrate its ability to perform adaptive mesh refinement using the DPG method's built-in error estimator. The final part of the study showcases the capabilities of the multiphysics framework through simulations of microwave-heated flow around obstacles with singular geometric features. These results highlight the potential of the proposed framework for the simulation and optimization of microwave-assisted chemical processes. Finally, the developed high-order multiphysics solver has a low memory footprint, since the electromagnetic solver relies on a Conjugate Gradient (CG) iterative solver and the fluid solver is implemented in a matrix-free fashion, making the overall approach scalable and well-suited for large-scale parallel simulations.

Nonlinear-manifold reduced-order models for parametrized finite element problems can achieve substantial compression both in the number of generalized (latent) coordinates and, through sampling-and-weighting hyperreduction, in the number of sampled elements/integration points. Yet current sampling-and-weighting approaches employ weights that remain fixed over the solution manifold. We contend that this restriction leaves hyperreduction potential untapped: allowing the weights to vary continuously and nonlinearly with the latent coordinates can further decrease the number of sampled spatial entities. To exploit this possibility, we propose the Manifold-Adaptive-Weight Empirical Cubature Method (MAW-ECM). Starting from a feasible fixed-weight ECM rule, a greedy pruning strategy removes sampled entities through convex quadratic weight-redistribution problems enforcing local conditions and positivity. The method is assessed on two nonlinear benchmarks: homogenization of a metamaterial unit cell exhibiting negative incremental stiffness, and a history-dependent continuum-damage problem. In both cases, the nonlinear manifold is constructed from an initial linear compression followed by an input-informed identification of the latent coordinates as general linear combinations of the retained modal amplitudes, incorporating graph information when relevant to seek the intrinsic dimensionality of the solution manifold. We show that combining the nonlinear-manifold representation with MAW-ECM reduces the number of sampled integration points by more than two orders of magnitude relative to the corresponding standard linear reduced model. Relative to the fixed-weight manifold models alone, the adaptive weights eliminate approximately 80% of the remaining points in the homogenization benchmark and more than 97% in the damage benchmark, while essentially preserving their accuracy.

Learned graph simulators provide an efficient alternative to high-fidelity solvers for granular dynamics. However, granular motion depends strongly on inter-granular contact history, which is difficult to preserve when particle contacts form, break, and rearrange. Existing simulators mainly store temporal information in node features or node-level memory. Here we introduce TRACE, a graph-network simulator that stores interaction history directly on contact edges. Each edge maintains a persistent memory updated by attention-based message passing and a gated recurrent unit, while an edge-identity dictionary preserves this memory as the contact graph changes. A physics-structured decoder predicts inter-granular normal and tangential contact forces, enforces the Coulomb friction limit, and applies equal-and-opposite internal forces. The model is trained with single-step pretraining followed by autoregressive rollout fine-tuning. We evaluate TRACE on 2D and 3D granular column-collapse benchmarks. In both cases, TRACE produces stable, physically consistent long-horizon rollouts, closely reproducing the final deposit geometry and the kinetic energy released during collapse. Compared with graph network simulator (GNS) and node-memory graph neural simulator (NMGNS), TRACE reduces long-rollout position error by 31-62% and final-deposit error by 58-89% across the two benchmarks, while using fewer parameters and maintaining near-zero particle interpenetration. TRACE also achieves 12.2$\times$ and 8.9$\times$ speedups over the material point method (MPM) reference solver in 2D and 3D, respectively. Our code is available at https://github.com/Data-Driven-Computational-Geotechnics/TRACE.

We prove that fully independent SparseStack achieves the oblivious subspace embedding parameters conjectured by Nelson and Nguyen (FOCS 2013): $m=O((d+\log(1/\delta))/\varepsilon^2)$ rows and $s=O(\log(d/\delta)/\varepsilon)$ nonzero entries per column for distortion $\varepsilon$ and failure probability $\delta$ on any fixed $d$-dimensional subspace, with explicit constants. The proof turns random-matrix concentration into a problem in finite-dimensional linear algebra. A coupling first reduces the moment estimates to a model with independent finite-valued entries. We represent these variables by multiplication operators, so their matrix moments become exact matrix elements of a deterministic operator on a finite tensor product. The central estimate bounds the contribution of $\ell\ge1$ occupied sites sharing the external factor $\mathbb{R}^d$ by $d+\ell-1$ rather than $d\ell$, yielding additive dependence on the dimension and the moment order. This approach controls both spectral edges without Gaussian comparison. The main theorem has been formally verified in Lean 4.

We develop an efficient and accurate multiscale computational framework for predicting the buckling and post-buckling behaviour of elastic beams containing periodically distributed soft inclusions. The framework extends our previous multiscale, patch, computational homogenisation for linear elasticity by incorporating nonlinearity, and thus enables accurate prediction of both the buckling onset and the subsequent post-buckling response. Microscale computations are performed only within a sparse set of small subdomains (patches), while the macroscale behaviour is recovered through a proven patch-coupling algorithm. The scheme is assessed through quarter- and half-domain patch computations for beams with inclusion-to-matrix Young's modulus ratios ranging from 0.001 to 1. The results show that reducing the inclusion stiffness lowers both the critical buckling strain and the critical buckling stress, indicating an increased susceptibility to instability, while producing a milder post-buckling response with smaller transverse deflections and stress drops. Eigenvalue analysis of the Jacobian matrix accurately predicts the onset of instability and the corresponding critical strain and stress. Bifurcation diagrams of the nonlinear buckled configurations under compressive loading, and a quantitative analysis of the effect of the interpolation order on the predicted buckling and post-buckling responses, are also presented. Comparisons with full-domain simulations demonstrate that the proposed framework accurately predicts both the buckling threshold and the post-buckling behaviour while substantially reducing the computational cost. The methodology is readily extendable to heterogeneous beams, plates, shells, and other engineering structures.

Ordinary differential equations (ODEs) underlie models in science and engineering, and many applications need derivatives of their solutions with respect to parameters. Ensembles of independent trajectories suit graphics processing units (GPUs), but current GPU software forces a trade-off: the fastest ensemble solvers cannot be differentiated in reverse mode at the speed they solve, and the solvers built for differentiation solve more slowly. No single tool has yet offered a reverse-mode gradient at the speed of a fused-kernel solve. We present GRADSOLVE, an open-source JAX library for solving and reverse-mode differentiating low-dimensional ODE ensembles on NVIDIA GPUs. It records the steps an adaptive solver accepts and differentiates a fixed-step replay of them; the returned gradient is the exact discrete adjoint of those steps, the same derivative Diffrax returns by default, obtained more cheaply from a fixed-length chain than from an adaptive loop. It targets ensembles differentiated many times against one recorded mesh, keeps Diffrax as a fallback, and supports explicit and Rosenbrock integrators. Used as a solver, GRADSOLVE's forward-only kernel ran 2.8x faster than DiffEqGPU.jl; used for gradients, once a record exists, it computed them 5.6-14.1x faster than Diffrax's checkpointed adjoint at matched forward-state accuracy across three GPU generations, the advantage narrowing on large ensembles and, on stiff systems, down to parity at tight accuracy. GRADSOLVE is released at https://github.com/ECLIPSE-AI4Science/gradsolve.

Recent advances in equation discovery methods such as SINDy have highlighted the growing interest in identifying governing parameters and models directly from data. In this work, we take a complementary approach grounded in analysis and numerical PDE methods: we recover an unknown viscosity in steady Navier-Stokes equations (NSE) from partial incompressible flow observations using continuous data assimilation (CDA). We propose a simple and efficient parameter recovery algorithm and also a nonlinear solver for CDA-NSE. Together, this creates a highly efficient technique for recovering an unknown viscosity from partial solution data. Our analysis establishes the well-posedness of steady CDA-NSE, quadratic convergence of the parameter recovery algorithm, and quadratic convergence of a CDA-Picard + CDA-Newton nonlinear solver. Numerical experiments illustrate that the methods are very effective in restoring parameters quickly, even with poor initial guesses.

We approximate strongly continuous convex monotone semigroups by learning their Chernoff-type one-step operators with neural operators. First, we introduce the general class of so-called Chernoff-neural operators and show in a universal approximation theorem that they can approximate the Chernoff one-step operators arbitrarily well. By using stability estimates between weighted Hölder spaces, the one-step approximation error can be propagated through the iterations which yields universal approximation of the corresponding semigroup. Second, we introduce the more specialized class of envelope-neural operators for envelope semigroups which allows us to derive quantitative approximation rates. Finally, we illustrate the effectiveness of these neural operators in several numerical examples arising from non-linear partial differential equations, stochastic optimal control and stochastic processes under model uncertainty.

We consider the numerical solution of the linear Vlasov transport equation on bounded spatial domains with inflow boundary conditions based on low-rank approximation. We combine a finite volume discretization with a rank-adaptive step-and-truncate scheme for the resulting matrix ODE. The spatial and velocity meshes may be unstructured, while suitable numerical fluxes retain a separated space-velocity representation. For homogeneous inflow, we show that the low-rank scheme inherits the $L_2$ stability and CFL restriction of the underlying full finite volume forward Euler method. In addition, the low-rank approximation error is bounded explicitly in terms of the truncation tolerances, avoiding the modeling error associated with tangent-space projections in dynamical low-rank approximation. Numerical experiments in 1d1v and 2d2v confirm the predicted error behavior. In 2d2v, the method handles an unstructured spatial mesh with nonzero inflow and a full tensor-product discretization of approximately $5.8\cdot10^{10}$ phase-space cells while the numerical rank is at most twelve.

How small can the positivity loss be for an $N$-bandlimited spectral operator that exactly reproduces all modes up to degree $L$? For spherical polynomial approximation on $\mathbb S^d$, we prove that the smallest possible excess of the uniform operator norm above 1, equivalently the least positivity loss, is of sharp order $\left({L}/{(N+1)}\right)^2$ when $1\le L< N$. The lower bound follows from a Fejér peak test and a concentration estimate for bandlimited kernels, while a matching upper bound is obtained by correcting a positive Jackson operator with a smooth filter. We illustrate the result in three settings. On the circle, taking $N=sL-1$, this determines the sharp order of the generalized-projection constant above 1 and identifies the gap between the $s^{-1}$ excess of delayed de la Vallée–Poussin means and the optimal $s^{-2}$ order. For filtered hyperinterpolation, whose operator norm has long been known to be uniformly bounded, we give a quantitative lower bound on its separation from the positivity threshold 1. Finally, we identify an operator-level obstruction to maximum principles.

This paper studies the computational implementation of proper orthogonal decomposition reduced-order models (POD-ROMs) for evolutionary parametric time-dependent partial differential equations (PDEs). For a one-parameter model, many papers in the literature build the correlation matrix by projecting onto $L^2(\Omega)$ even though optimal pointwise error estimates are proved when projecting onto $H^1_0(\Omega)$. We compare both scenarios and observe the similar performance in practice. Additionally, to get the POD approximation it is necessary to compute the nonlinear term in the reduced equations. This requires a high computational cost that increases when the problem gets more complex. Numerical results in this paper show different approaches to address this issue. Discrete Empirical Interpolation Method (DEIM) is the most efficient approach. It reduces computational time by approximately half compared to computing the whole FEM formulation by a tensor construction while maintaining the same accuracy. For a multiparameter model, we analyze the new and standard method proposed in \cite{newmethod} using a two-dimensional Brusselator model to complement the results and support the error analysis with a more complex system.

The Boltzmann equation plays an important role in modeling mesoscopic behavior in a wide range of scientific and engineering applications. However, its numerical solution is computationally expensive due to the high dimensionality of the model and the nonlinear nonlocal collision operator, especially for steady-state problems that require iterative solvers. This cost becomes prohibitive for inverse problems, where the induced optimization problem requires repeated forward solves. In this work, we propose a reduced-order model (ROM) for the parametric Boltzmann equation to address this computational challenge. The ROM constructs a low-dimensional approximation space for the parameter-induced solution manifold through a residual-based greedy strategy, and the reduced solution is then obtained via residual minimization over the reduced space, subject to mass conservation. The overall efficiency of the ROM is achieved by exploiting the quadratic structure of the collision operator and a precomputed separable approximation of the collision kernel. The resulting ROM is further applied to a thermally-driven inverse problem for reconstructing collision parameters from the observed macroscopic temperature data. This is accomplished either by directly replacing the PDE constraint with the ROM, leading to a bilevel optimization formulation, or by reformulating the task as a single-level optimization problem through the Karush–Kuhn–Tucker (KKT) conditions. Numerical experiments in both collision-dominated and transport-dominated cases are performed to demonstrate the efficiency and accuracy of the proposed ROM and its effectiveness in inverse problems. In particular, the resulting inverse problem is computationally much more tractable, achieving speedups of several orders of magnitude over that based on the full-order model while maintaining comparable accuracy.

In many systems, the interaction between a deformable surface or interface and the surrounding bulk fluid is coupled with intrinsic spatiotemporal dynamics within the moving surface. Examples include tumor growth, biological tissue morphogenesis, cardiac mechanics, multi-phase surfactant chemistry, additive manufacturing, clothing wear-and-tear, and reactive combustion flows. Solving such problems requires both geometric computing algorithms to track and resolve the surface and numerical methods to solve the coupled governing equations in the surface and the surrounding bulk phase. Here, we present a fully meshfree numerical solver for such coupled bulk-surface problems with deformable interfaces. The presented solver tracks the surface implicitly, solving for the dynamic surface geometry based on stress balance coupled to surrounding fluid phases. We show convergence for a mass-conserving case on a growing sphere and solve bulk-surface problems with incompressible Navier-Stokes fluids coupled to in-surface nonlinear reaction-diffusion dynamics. Finally, we show a model of biological morphogenesis, solving simultaneously for the dynamic surface shape and the fields on the curved surface with two-way coupling.

In this paper, we consider the Smagorinsky model for the Navier-Stokes equations within a virtual element framework. Under the standard assumption of small data, we prove the existence and uniqueness of a solution. Assuming more regularity to the solutions, we derive the known convergence rates $h$ for the a priori error estimates of the Smagorinsky model in two dimensional domains. We additionally prove that divergence-free virtual discretizations provide improved convergence orders, with weaker regularity assumptions than in the finite element literature. We conclude the paper with numerical results that corroborate the theory.

Determining the dimension of the $C^1$ cubic spline space $S_3^1(\mathcal{T})$ on an arbitrary nondegenerate planar triangulation has remained unresolved since the 1970s. Schumaker's lower bound includes a local correction $\sigma$ for singular interior four-stars, and it was conjectured that this bound is always attained. We disprove this conjecture by constructing a one-parameter family of nondegenerate realizations of a fixed 18-triangle complex, with only the central vertex moving as $v_6(t)=(t,0)$ on the admissible interval $I=(-3/4,24/55)$. The family exhibits three distinct cases. For $t\in I\setminus\{1/5,3/83\}$, the lower bound is attained and $\dim S_3^1(\mathcal{T}(t))=33$. At $t=3/83$, the central four-star is singular, $\sigma=1$, and the resulting dimension 34 is exactly accounted for by the classical local correction. At $t=1/5$, however, all interior vertices are nonsingular and $\sigma=0$, yet $\dim S_3^1(\mathcal{T}(1/5))=34>P_{\mathcal{T}(1/5)}(1,3)=33$. The smoothing-cofactor calculation shows that the dependence at $t=3/83$ is confined to the central vertex block, whereas the dependence at $t=1/5$ couples all seven interior vertex cycles even though every individual block has full row rank. A complementary Bernstein–Bézier calculation gives the same dimension profile. Thus the singular-four-star correction does not capture every geometry-dependent contribution to $\dim S_3^1(\mathcal{T})$; genuinely global compatibility must also be taken into account.

We consider reaction-diffusion eigenproblems with oscillatory diffusion and reaction coefficients. The reaction coefficient magnitude is large: the corrector equation identified by periodic homogenization involves both the diffusion and the reaction operators. We study the numerical approximation of this problem using the Multiscale Finite Element Method (MsFEM). This now classical method is a finite element type method that performs a Galerkin approximation of the oscillatory problem on a specific, problem dependent, basis set. The basis functions are precomputed in an offline stage. Inspired by homogenization theory and using some filtering ideas, we show how to define these basis functions in order to obtain an efficient method. The comprehensive set of numerical experiments that we present, in periodic and non-periodic cases, for the scalar-valued version of the problem (which is then self-adjoint) and for the vector-valued version of the problem (which is then in general non self-adjoint), demonstrates the performance of the approach. Some theoretical arguments complement the numerical observations.

We introduce a stream-function formulation for the Surface Stokes-Helfrich model and study an Arbitrary Lagrangian-Eulerian (ALE) Surface Finite Element (SFEM) discretization in space and a semi-implicit discretization in time of these equations. In contrast with special cases of a surface Stokes model on stationary or prescribed evolving surfaces the surface pressure remains as an unknown and needs to be reconstructed in each time step. As for the surface Stokes model on a prescribed evolving surface the formulation also requires an additional gradient potential. A detailed computational comparison with an ALE-SFEM discretization in space and a similar semi-implicit discretization in time of the corresponding velocity-pressure formulation is explored on a representative problem formulation for simply-connected surfaces, demonstrating the validity of the approach.

We introduce a compactly supported multinode Shepard operator for the interpolation of scattered data on the torus embedded in $\mathbb{R}^3$. The method combines local polynomial interpolation of total degree $d\in \mathbb{N}$ with compactly supported Shepard-type weights, so that the approximation at each evaluation point depends only on neighbouring stencils of nodes. The torus is treated as an algebraic surface defined by a quartic polynomial, and Gröbner bases are used to construct reduced polynomial spaces on the surface by removing the redundancy induced by the defining equation. This yields local Vandermonde systems adapted to the toroidal geometry. We discuss the metric structure of the torus, show the local equivalence between the periodic parameter distance and the Euclidean distance inherited from $\mathbb{R}^3$, and use this equivalence to motivate the compact support construction. We establish a uniform error estimate in terms of the maximal support radius and the local Lebesgue constants. Under uniform locality and stability assumptions, the method converges with order $d+1$ with respect to the fill distance. Numerical experiments on analytical test functions confirm polynomial reproduction and exhibit an error decay consistent with the theoretical analysis. The approach is also tested on Computational Fluid Dynamics data mapped onto the torus, including the interpolation of the velocity components, the reconstruction of the velocity magnitude from the interpolated components, and the reconstruction of a tangent velocity field through an orthonormal lifting of the interpolated components.

Numerical experiments reported by Quan and Wu [SIAM J Numer Anal 61 (2023) 2106-2132] show that the observed temporal convergence rates of nonuniform L2 methods for subdiffusion models are not consistent with the theoretically predicted order $3-\alpha$. This discrepancy suggests that a more refined analysis is needed and motivates the development of a nonuniform BDF3-L2 method for the subdiffusion equation. To account for the initial solution singularity, we employ the multiscale solution decomposition to decompose the original solution and approximate a smoother unknown variable that satisfies the subdiffusion model with a smoother source term. The resulting formulation, however, involves restrictive high-order boundary conditions on the source term and initial data. To overcome this difficulty, we introduce a spectral truncation technique that requires only slightly stronger regularity of the data and a controllable truncation error. We establish high-order regularity estimates of the solution to the truncated problem and develop a nonuniform BDF3-L2 method for its numerical approximation, based on which we derive a rigorous error estimate of temporal convergence order $2+\alpha$. Numerical experiments are carried out to substantiate the theoretical findings.

The Johnson-Lindenstrauss (JL) lemma guarantees that a random projection of $n$ points to $m=O(\varepsilon^{-2}\log n)$ dimensions preserves pairwise squared distances within relative error $\varepsilon$ with high probability, and this dimension order is asymptotically optimal. In high dimensions, however, distances concentrate around a baseline while key geometric information lies in much smaller fluctuations. We show that the JL bound can therefore be uninformative about retained geometry: an independent Gaussian replacement map can satisfy it even though the replacement cloud is independent of the original data. We then ask how well any decoder can recover a feature $f(D)$ of a squared distance $D$ from a linear sketch. Under squared-error loss, the optimal decoder is conditional expectation, so recovery defines a linear operator whose singular values quantify feature recovery. For isotropic Gaussian data ($\Sigma=\sigma^2 I_d$), we diagonalize this operator in closed form. For fixed $k$ with $m,d-m\to\infty$, its $k$th singular value satisfies $\ell_k\approx(m/ d)^{k/2}$. This yields three sharp consequences. A rank-$m$ sketch retains at most an $m/d$ fraction of the variance of any feature of one squared distance. If $m\to\infty$ and $m/d\to0$, the expected Kendall correlation is $\frac{2}{\pi}\sqrt{m/d}(1+o(1))$; for fixed $q$, nearest- neighbor agreement tends to $1/q$. Yet one projection can satisfy the JL bound while mean Kendall correlation vanishes when $\log n\ll m\ll d$. After removing scale, Haar-averaged retained covariance-shape information is $(m/d)^2$. Thus JL distance preservation does not quantify the geometry available for comparison or inference.

The Fokker-Planck equation admits both a macroscopic Eulerian description through probability densities and a microscopic Lagrangian description through stochastic trajectories. Consequently, optimization problems constrained by the Fokker-Planck equation can be formulated from either perspective. Surprisingly, the corresponding adjoint equations appear to be fundamentally different: the macroscopic adjoint is governed by the backward Kolmogorov equation, whereas the microscopic adjoint evolves pathwise along stochastic trajectories. In this note, we reconcile these two formulations by establishing their correspondence in the continuum setting. We further show that, although their discrete gradients no longer coincide after discretization, both provide consistent numerical approximations of the continuum gradient. Explicit convergence rates are established for both discretization strategies.

Kolmogorov-Arnold Networks (KANs) place learnable B-spline activations on network edges rather than fixed activations on nodes. The standard Cox-de Boor recursion evaluates these activations through $k$ sequential passes for degree-$k$ splines, consuming over 90% of forward-pass time. InKAN replaces this recursion with the truncated power form, a classical result from approximation theory that expresses each uniform cubic B-spline as five $(x)_+^3$ terms at shifted knot positions. The resulting expression computes exact B-spline basis values: the same mathematical function as the Cox-de Boor recursion, evaluated without sequential passes. This paper documents three contributions: (1) an implementation structured for torch . compile fusion, eliminating recursion, span lookup, and scatter-gather operations; (2) a bounded-coordinate evaluation that clamps the normalized input to $[0, k{+}1]$, preventing the growth of cancellation error at large off-support coordinates; and (3) an open-source package (pip install inkan). In the tested configurations, InKAN has 2.8–3.5$\times$ lower forward-pass latency than the Cox-de Boor recursion. Partition-of-unity errors remain below $10^{-5}$ for grid sizes up to 200.

The Jacobi-Davidson method is a widely used subspace method for computing a few eigenpairs of a large, sparse, non-Hermitian matrix closest to a target. Like other subspace methods, it orthogonalizes each new expansion vector against the whole search basis, at a cost that grows quadratically with the subspace dimension and, in a distributed setting, requires a global synchronization at every iteration. We introduce a randomized Jacobi-Davidson method that replaces this orthogonalization with a much cheaper randomized orthogonalization process, which requires no inner products involving full-dimensional vectors. We prove that, under a uniform separation condition, the randomized method retains the local quadratic convergence of classical Jacobi-Davidson for non-Hermitian problems. The key step is showing that an exact solution of the randomized correction equation is one step of sketched inverse iteration for an arbitrary shift, which lets us extend the analysis to harmonic and refined variants of the sketched extraction, better suited to interior eigenvalues. Numerical experiments on real and synthetic non-Hermitian eigenvalue problems confirm the predicted convergence order and show that the randomized method matches the reliability of classical Jacobi-Davidson while reducing orthogonalization cost.

Finite elements provide accurate and efficient methods for the numerical solution of partial differential equations by means of restricting variational problems to finite-dimensional approximating spaces. However, they do not, in general, guarantee enforcement of bounds constraints inherent in the original problem. We propose two approaches to enforcing bounds constraints for time-dependent problems. First, we propose general projective methods which result from a systematic modification of any abstract time-stepping scheme. Second, we present a monolithic technique for which we take a modified formulation of implicit single-stage Runge-Kutta methods and general implicit multistep methods as prototypical examples. By solving a constrained optimization problem, we are able to ensure that the bounds constraints are enforced by the approximate solution at the discrete time levels, obtaining (formally) high order methods both in space and time. Numerical examples for the linear heat and advection equations and nonlinear Allen-Cahn equation are given.

We develop a structure-preserving, data-driven collision operator for spatially homogeneous multi-species kinetic systems from molecular dynamics (MD). The operator consists of diagonal self-collision blocks and ordered off-diagonal cross-species blocks to describe intra- and inter-species momentum and energy exchange. Within a local and point-wise identifiable kernel class, we develop the necessary and sufficient condition for the admissible kernel class satisfying the conservation laws, the H-theorem, and the frame indifference. Unlike the classical Landau operator, the off-diagonal kernels are not restricted to be symmetric under permutation of the two velocity variables. This unique structural freedom captures the distinct responses of different species to unresolved correlations and many-body effects arising from micro-scale particle interactions. The equivalent parameterizable kernel formalization enables us to learn a generalized data-driven collision operator directly from MD, where the low-rank tensor representations and random sampling are used to achieve efficient kernel training and numerical simulation. Numerical experiments show that the learned operator accurately predicts transport coefficients and the non-equilibrium relaxation, while retaining discrete conservation and entropy production. In particular, it captures plasma kinetics in the moderately coupled regime, where the predictions of both the Landau and the data-driven model restricted to velocity-permutation symmetry show significant discrepancies.

We study minimal residual space-time finite element discretizations of linear parabolic initial value problems in canonical space-time variational form. To deal with the arising dual norm, we introduce the Riesz lift of the residual as an additional variable. Quasi-optimality of the primal variable of the mixed system follows from a uniform inf-sup condition. This condition is known to be satisfied for finite element spaces w.r.t. prismatic partitions of the space-time cylinder that allow for a decomposition into time-slabs. We prove that this condition cannot be expected to hold otherwise. To recover stability for general partitions and the data at hand, we derive an a posteriori condition on the error between the exact Riesz lift of the residual and its Galerkin approximation – being the secondary variable of our system – under which the primal variable is quasi-optimal. We derive a posteriori error estimators for both variables, and use them in a double-adaptive loop that alternates test-space with trial-space enrichment. We illustrate our findings with numerical experiments in $1+1$ and $2+1$ dimensions.