Papers

In this paper, we unravel the effect of different decoder architectures on the interpretability of the latent space of the Physically Recurrent Neural Network. Particular emphasis is given to a new weight normalization constraint, which acts as a regularization technique and enables robust training in the low-data regime. A brief visual exploration illustrates how these changes impact the latent space and how the fictitious stress can align with the true state of the RVE without explicit training. Reaping the benefits of a meaningful latent space, a case study illustrates how information from the microscopic level can be retrieved and incorporated into a multi-task approach that does not require extra parameters or larger training sets. Another key contribution shows that a specific architectural choice can naturally lead to a thermodynamically consistent formulation. By enforcing an adjoint encoder-decoder structure with positive scalar contributions, this modification ensures energy consistency across scales and non-negative dissipation, leading to even lower training requirements. This alternative completes the study on interpretability, inductive bias, and thermodynamic consistency, and demonstrates that data efficiency can be improved with careful architectural choices rooted in the underlying physics.

We propose a data-driven framework for estimating the Bayesian Cramér-Rao bound (CRB) in high-dimensional imaging systems with complex, analytically intractable priors. Direct CRB computation is challenging in this setting due to the need to model the prior score and to form and invert the Bayesian Fisher information matrix in very high dimensions. To address these issues, we first reformulate the inverse problem in the latent space of a pre-trained variational autoencoder, thereby dramatically reducing the dimensionality of the bound estimation problem while preserving the spatial structure of the images. The Bayesian CRB is formed in this latent space and mapped back to the native parameter space using a change-of-variables formula. Second, to learn the latent prior score, we introduce a new sliced score matching objective defined in a Bochner space endowed with an $H^1(\Omega)$ Sobolev norm in space. This "Bochner-space sliced score matching" objective is consistent with standard sliced score matching, but penalizes errors in the spatial gradients of the score, suppressing high-frequency artifacts that otherwise contaminate the resulting CRB estimates. We validate the approach on a stylized quantitative photoacoustic computed tomography (qPACT) breast imaging problem with over one million unknown parameters, using a foundation-model autoencoder derived from Stable Diffusion. The proposed method yields stable, artifact-free Bayesian CRB estimates that reflect the highly non-Gaussian structure of the learned prior and reveal the substantial impact of the prior on the relative performance of competing qPACT design schemes.

We develop a Bayesian framework for learning nonlinear operators between infinite-dimensional spaces. Given a map $g_0:\mathcal{X}\to\mathcal{Y}$ between separable Hilbert spaces, we study the recovery of $g_0$ from $n\in\mathbb{N}$ noisy input-output pairs $(\boldsymbol{X},\boldsymbol{Z})=(X_i,Z_i)_{i=1}^n$ with $Z_i= g_0 (X_i ) + E_i$. Here the $X_i\in\mathcal{X}$ are randomly drawn 'design' points in a compact subset of $\mathcal X$, and the $E_i$ are assumed to be i.i.d. draws from a Gaussian white noise process indexed by $\mathcal{Y}$. For any 'operator-valued' prior $\mathbb{P}_G$ supported on the space of continuous operators, we show existence of the posterior $\mathbb P_{G|(\boldsymbol{X},\boldsymbol{Z})}$ as a regular conditional distribution, and provide a characterization of its Radon-Nikodym derivative. For Gaussian priors, we establish algebraic (in the sample size $n$) convergence rates for the posterior mean towards the ground truth; this corresponds to a ridge regularized kernel estimator. Moreover, we show that the posterior mean is minimax optimal (up to logarithmic factors) over hyperrectangles when the smoothness of the prior matches that of the ground truth. To illustrate the applicability of our analysis, we derive explicit learning rates for the Darcy flow solution operator.

We study a physics-informed neural network (PINN) approximation of a coupled PDE system that models the interaction between fluid and elastic dynamics. The canonical heat-wave model considered here couples a parabolic and a hyperbolic equation through transmission conditions on a fixed boundary interface, and serves as a prototype for the Stokes-elastic systems that arises in real world applications. This coupling makes the error analysis more delicate than that for a single PDE. We bound the error between the exact and PINN solutions by the associated PINN residuals, and we moreover prove that the generalization error is controlled by the training error, the number of quadrature points, and the size of the network. {To the best of our knowledge}, these are the first error estimates for a PINN approximation of a coupled fluid-structure interaction PDE system. We also show that numerical experiments with respect to a given test problem.

We extend to the Willmore flow the phase-field neural operator approach proposed by Bretin, Denis, Masnou, and Terii (2022) for mean curvature motion in two and three dimensions. The resulting lightweight, single-time-step architecture is motivated by a Bence-Merriman-Osher-type expansion of the Allen-Cahn semigroup, and is trained using trajectories generated by a hybrid minimizing movement scheme. In our experiments with non-oriented interfaces, unconstrained convolution kernels produce anisotropic dynamics, whereas radial kernels enhance the evolution process. We also construct three-dimensional kernels from trained two-dimensional kernels via their radial Fourier profiles, which provides an effective initialization for the non-oriented 3D model. Finally, we show applications to curve and surface reconstruction from unoriented point clouds, including reconstruction using the transferred model without the need for additional 3D training.

This paper develops a deep learning framework for solving high-dimensional stationary Fokker-Planck equations by combining Monte Carlo weak supervision, energy-variational physics constraints, and quasi-Newton refinement. Monte Carlo statistics anchor the scale of the stationary density, while a free-energy variational loss enforces equilibrium physics; the formulation also applies to stationary systems with pointwise orthogonal rotational components and, using a family of test functions, to systems with more general non-gradient drift. Numerical results demonstrate accurate and robust density approximations across a range of high-dimensional stationary Fokker-Planck problems, including rough-potential and noisy-data settings.

We propose and analyse a fully mixed finite element method for a two-way coupled mechanochemical model of calcium signalling in viscoelastic tissue. The mechanical response follows a Kelvin–Voigt law, with the elastic and viscous stresses, velocity, and spin tensor as unknowns, and the symmetry of the total stress imposed weakly. The reaction-diffusion system for calcium and inactivated inositol triphosphate is also formulated in mixed form, introducing the diffusive flux as an additional unknown. The coupling is bidirectional: calcium generates an active stress, while local dilation feeds back into the calcium dynamics. We establish global existence and uniqueness by reducing the coupled system to an integro-differential problem and subsequently recovering the mixed variables. The semidiscrete scheme is formulated for general conforming finite element spaces and shown to be well posed under a single inf-sup condition. An interpolation operator that commutes with the divergence and preserves the weak symmetry pairing yields optimal a priori error estimates. We specialise the analysis to three families of weakly symmetric finite elements, deriving explicit convergence rates and, for one family, a superconvergent local post-processing of the displacement. Numerical tests in two and three dimensions confirm the predicted convergence rates and illustrate the method through viscoelastic benchmarks and mechanically driven calcium wavetrains.

We present a space-time finite element formulation for the transient analysis of geometrically exact shear-deformable beams. The formulation treats space and time within a unified finite element framework and employs a mixed approach in which translational and angular velocities are introduced as independent fields. The rotation field is parametrized using nodal rotation (pseudo-)vectors and a suitable objective interpolation of the triad field is employed. The resulting formulation supports generally curved beam structures, including geometric kinks and rigid joints. Furthermore, time-dependent material domains can be represented directly through the geometry of the space-time surface, enabling the treatment of non-material boundaries without modifying the underlying beam formulation. Numerical examples demonstrate convergence to solutions obtained with classical time-stepping schemes, objectivity of the formulation, and the influence of the space-time stabilization. Finally, the sliding spaghetti problem illustrates the applicability of the proposed approach to axially moving continua.

This article examines the Fictitious-Domain approach using Distributed Lagrange Multipliers (FD-DLM), a method that handles elliptic interface problems with separate meshes without requiring the interface to match the background mesh. It compares two stable, low-order FD-DLM discretizations when using both continuous and discontinuous multiplier spaces. The study distinguishes between satisfying the weak coupling constraint per element and the size of the discrepancy between the two discrete primal fields. Numerical experiments indicate that the continuous formulation usually yields smaller local errors and maximum mismatch values, while the discontinuous formulation shows a more regular decrease in mismatch as the mesh is refined and has smaller primal errors.

We construct an explicit quartic polynomial for which Newton's and Halley's root-finding methods both fail to converge to a root on nonempty open sets of initial values. The example can be found through a reproducible computational search within a one-parameter family of quartics designed to exhibit a prescribed Newton two-cycle. We then prove analytically that Newton's method has a superattracting two-cycle and that Halley's method has a distinct attracting extraneous two-cycle, using a two-interval contraction argument. This fills a degree-four gap in the simultaneous Newton–Halley failure problem. The computational approach also identifies another quartic with numerical evidence of simultaneous failure of Newton's, Halley's, and Schröder's methods.

In this paper, we study the problem of sampling from a probability distribution using Langevin-type dynamics. We consider overdamped Langevin dynamics perturbed by adding a Stratonovich perturbation that preserves reversibility. We prove an optimal scaling for the Stratonovich perturbation applied to Gaussian target distributions and derive an algorithm to construct such perturbation. Our theoretical results are supplemented by numerical experiments in which we compare our proposed method with nonreversible samplers and overdamped Langevin dynamics. Our numerical experiments demonstrate the efficiency of the approach beyond the quadratic case.

Positive reciprocal pairwise comparison matrices are studied in entrywise logarithmic coordinates. Starting from the established consistent-residual split, we show that, for any fixed ordering, the residual matrix admits a unique top-down decomposition into mutually Frobenius-orthogonal layers. The construction applies to every positive reciprocal PC matrix of size $n \geq 3$, without additional spectral, sparsity, or genericity assumptions. An explicit active-matrix recursion and an equivalent cumulative parameter-extraction procedure have computational complexities $O(n^3)$ and $O(n^2)$, respectively. The squared Frobenius norms of the layers define an ordering-conditioned residual-energy profile that provides structural information complementary to the geometric-mean priorities and the Koczkodaj inconsistency indicator (Kii). Entrywise exponentiation yields an exact Hadamard-product reconstruction of the original reciprocal matrix. Numerical experiments confirm agreement between the two computational procedures, exact reconstruction, mutual orthogonality, and numerical accuracy to machine precision.

We develop and analyze a physics-driven multiscale method for model reduction of the Helmholtz equation at high frequency. The method constructs a reduced approximation space from local wave responses within a mixed finite element formulation with velocity elimination. These responses are computed by solving Helmholtz problems with boundary excitations on oversampled regions, giving snapshots that incorporate the oscillatory behavior and propagation characteristics of the underlying medium. Local generalized eigenvalue problems then select modes from this physically informed space to form a reduced global system with sparse interelement coupling. To establish stability and convergence of this reduction, we relate the local spectral approximation to approximation of the global adjoint problem. Starting from an inf–sup estimate for the fine-grid formulation, this argument yields a multiscale inf–sup condition governed by fine-grid resolution, oversampling, and the first omitted local eigenvalues. The resulting stability establishes well-posedness and quasi-optimality of the multiscale approximation in a discrete energy norm. Together with local spectral estimates, it gives an a priori error bound that explicitly quantifies the roles of the coarse and fine mesh sizes, oversampling depth, spectral truncation, and wavenumber. This analysis provides a quantitative link between local space reduction and the stability and accuracy of the global approximation. Numerical experiments in two and three dimensions demonstrate accurate wavefield approximation with substantial reductions in global dimension and solution cost for homogeneous and strongly heterogeneous media, including perfectly matched layer truncations.

Direct approaches to discrete parabolic Carleman estimates often rely on weighted identities tailored to the discrete operator. We prove an abstract transfer theorem that derives such estimates for fully discrete backward Euler schemes with conforming spatial reconstructions from a continuous Carleman inequality together with suitable weighted approximation and compatibility properties. For a class of time-independent positive self-adjoint operators, the key construction assigns to each discrete state a shifted continuous auxiliary problem whose discrete solution is exactly the prescribed state. Weighted space-time error estimates, uniform over the full discrete load space, then transfer the inequality while preserving the residual and observation terms of the original scheme. This separates the PDE-specific Carleman analysis from the numerical verification required for each discretization. We verify the hypotheses for the heat equation discretized by consistent-mass $P_1$ finite elements on locally graded meshes in two and three dimensions and by standard Cartesian finite differences in any fixed dimension. For the heat equation, the transferred estimate is valid for Carleman parameters up to the scale $\min\{h^{-4/5},(\delta t)^{-2/5}\}$. For bounded real-valued space-time potentials, the resulting estimates yield relaxed observability and approximate null controls with uniformly bounded cost and exponentially small terminal errors under independent spatial and temporal refinement; weak accumulation points of the reconstructed controls are continuous null controls.

We present a root isolation algorithm for an analytic function. The method combines interval subdivision with a heuristic test based on local cubic Hermite interpolation. Whenever the function is sufficiently close to its cubic Hermite interpolant at a few sample points, we assume that the critical points of the cubic interpolant give a valid root isolation for the original function. Such an algorithm cannot, in general, provide a rigorous certification of the computed intervals. Nevertheless, extensive numerical experiments show that it is highly efficient and almost always produces the correct results. We hope that a probabilistic or exact algorithm can be derived from this work.

Consensus-based optimization (CBO) is a robust and flexible zero-order, multi-particle class of methods devised to provide provable solutions to high-dimensional global optimization problems, including truly nonconvex and/or nonsmooth cases. The algorithmic framework is grounded in a principled trade-off between stochastic exploration of the search space and deterministic contraction toward a consensus point. This thesis is concerned with three recently established variants of CBO. The $\delta$-CBO method is similar to the classic CBO method, with the difference that the diffusion function is constant. The consensus freezing (CF) scheme can be understood as a piecewise CBO method. The overall time horizon is divided into intervals in which the consensus point is fixed. The consensus hopping (CH) scheme can be understood as a deterministic iterative scheme that jumps at each iteration to a Gaussian distribution centered at the consensus point of the previous distribution. The central contribution of this thesis is the rigorous derivation of propagation of chaos results for all three schemes. For each variant, we establish that the empirical measure of the interacting particle system converges, as the number of particles $N$ approaches infinity, to the solution of the corresponding mean-field equation. Two regimes are treated: first, on a high-probability set on which the empirical measures are uniformly bounded; second, in the unrestricted setting, where moment bounds and a careful treatment of the inter-particle dependence yield convergence without imposing conditioning on favorable events. The propagation of chaos estimates are derived by the coupling method. Together, these results provide quantitative propagation of chaos estimates for all three CBO variants and establish that the particle systems faithfully approximate the idealized mean-field dynamics in the large-particle limit.

In this paper, we propose a highly efficient and accurate method, called “lift-and-decoupling generalized transferable neural network" (LD-GTransNet), for numerical solution of elliptic and moving interface problems. Our method makes use of the multi-layer GTransNet approach and incorporates a novel {lift-and-decoupling} strategy, in which the hidden-layer neuron parameters are preset and the physical coordinates and the added auxiliary variable are processed by separate subnetworks. Such design enables a unified global approximation of piecewise smooth solutions without explicit domain decomposition, while preserving network expressivity and avoiding the redundancy inherent in conventional {lift} approaches. The resulting formulation reduces to a linear least-squares problem, thereby eliminating the need of nonlinear training and ensuring the high computational efficiency and accuracy. Extensive numerical experiments on benchmark interface problems in two and three dimensions demonstrate that, compared with existing neural network-based and traditional numerical methods, our proposed LD-GTransNet consistently delivers superior performance, particularly for problems with high-contrast coefficients and complex interface geometries.

Inverse design has been transformative in nanophotonics, providing an automated means of realizing high-performance, non-intuitive devices. However, the primary bottleneck of existing inverse design approaches is their reliance on time-consuming and computationally expensive full-wave electromagnetic simulations. In this work, we extend the Precomputed Numerical Green Function (PNGF) method for the first time to nanophotonic inverse design, accelerating evaluation of the forward problem by multiple orders of magnitude with no loss in accuracy compared to the full-wave solution. After a single, fully parallelizable precomputation step, both the objective function and its gradient can be evaluated through simple calculations whose cost scales linearly with the size of the design region. A low-rank matrix update technique further reduces the cost of objective function evaluation, yielding sub-millisecond computation times per iteration. Design examples using both gradient-based level-set optimization and tile-flipping direct binary search often achieve more than three orders of magnitude speedup in design time, resulting in an ultrafast photonic inverse design platform that can design new devices from scratch in seconds.

We develop a method for lossy compression of incompressible (i.e., divergence-free) velocity fields that preserves zero divergence and, on periodic domains, periodicity. We fit the nodal velocities with property-preserving matrix-valued kernels and compress their coefficients using standard lossy compressors. Every decompressed coefficient vector defines a field with the same constraints, independently of the compression tolerance. We present several fitting schemes designed around either low-rank representation of or sparsity in the interpolation matrix. More specifically, we show results for global Matérn interpolation, a two-level Nyström approximation, Nyström–Wendland residual interpolation, multilevel Wendland interpolation and a property-preserving (redundant) kernel frame. We derive nodal error bounds that relate the compressor tolerance to the additional velocity error. Experiments with four compressors on cavity, Taylor–Green, and species-transport flows verify the bounds under their coefficient-error assumptions and compare storage at matched velocity accuracy. The low-rank methods give compact representations when they resolve the velocity field sufficiently accurately; residual, multilevel, and frame interpolation extend the attainable accuracy when they do not. We also measure changes in turbulent velocity statistics. The resulting representations support an accuracy–compression tradeoff while preserving incompressibility and periodicity analytically.

Inferring an Allen-Cahn constitutive force from images requires more than an accurate optimizer: nonlinear observation noise alters the mean of force features, shared noise couples instruments to the response, and a temporal balance borrowed from another solver can remain biased even on clean data. We develop observation-aware, noise-adjusted moment equations that address these effects before inversion. The surface-tension-calibrated constitutive transformer (SCCT) couples the resulting moments to a positive Bernstein inverse and can adapt its prior and precision from images and moments. A conditional stability bound separates moment residuals from prior error. Synthetic tests show that correcting the estimating equation matters more than increasing regularizer flexibility; the benefit of learned adaptation is smaller and depends on the training comparison. Forces inferred from planar movies predict unseen planar and three-dimensional evolutions. A ring-cage breakup shows why diffuse-interface accuracy distinguishes forces even when their topological predictions agree, while external surface tension sets the energy scale separately.

Magnetic resonance elastography (MRE) is a noninvasive imaging modality for quantifying the viscoelastic properties of soft tissues from shear wave propagation. Recovering the complex-valued shear modulus from measured displacement fields leads to a severely ill-posed inverse problem, particularly in the presence of noise and limited boundary excitations. We investigate adjoint-based optimization, convolutional neural network (CNN) reconstruction, and a hybrid framework combining both approaches. The forward model is based on a scalar form of the modified stationary Stokes system with a complex shear modulus. We establish well-posedness of the forward problem, existence of minimizers, and first-order optimality conditions for the adjoint-based formulation, and implement a nonlinear conjugate-gradient method with Armijo line search. While PDE-constrained optimization can accurately refine coefficient reconstructions, its performance depends strongly on initialization. We therefore construct a two-dimensional CNN that maps complex-valued displacement measurements to spatially varying complex shear modulus fields and provides rapid, informative initial reconstructions. The proposed hybrid method uses the CNN reconstruction to initialize the adjoint-based optimization, yielding faster convergence and improved accuracy. The CNN is trained on coefficient fields containing individual perturbations and tested on both individual and previously unseen combined configurations. Numerical experiments demonstrate that the CNN generalizes to these more challenging configurations, while subsequent PDE-constrained optimization further refines the reconstructed coefficient. These results demonstrate the potential of combining data-driven initialization with physics-based optimization for efficient and accurate MRE reconstruction.

We propose a two-level filtered subspace iterative (TLFSI) method for computing all eigenpairs whose eigenvalues in the interval of interest $[a,b]$ after finite element discretization of self-adjoint PDE eigenvalue problems. The method combines a two-level approach with a modern filtered subspace iteration framework. In this paper, we use contour integral filtering as an example. The two-level method constructs a coarse space according to the spectral location of the eigenvalues of interest and obtains from it an approximation to the target invariant subspace. We prove that the distance between the approximate and target subspaces depends only on the separation between the spectral region of interest and the remainder of the spectrum, together with the coarse mesh size \(H\). Numerical experiments on benchmark elasticity problems confirm the theoretical predictions.

We develop a closure modeling framework for identifying missing components of dynamical systems using Kernel Ridge Regression (KRR). The framework addresses two classes of closure problems: difference-equation closures arising in ODE and PDE settings, and algebraic closures arising from moment closure in kinetic equations. For the first class, we derive an error bound in an ODE setting that quantifies contributions from time integration, approximation of unresolved scales, and interpolation required to couple unresolved-scale effects to the resolved solver. Numerical experiments on the Lorenz-63 system and the Kuramoto-Sivashinsky equation demonstrate accurate long-horizon predictions and substantial improvements over an LSTM-based closure model. For the second class, we consider moment closure for a one-dimensional kinetic equation by modeling discrepancies between kinetic and macroscopic fluxes as a function of the resolved macroscopic variables. We compare global KRR models based on PCA coordinates with spatially local models. While the global model performs well for unimodal initial conditions, its accuracy deteriorates for bimodal initial conditions. Spatially local models with appropriate modeling inputs improve robustness and achieve higher predictive accuracy.

We study the universal isotypic block-diagonalization problem in the Type-2 theory of effectivity (TTE). We obtain a uniform computability result for the universal block diagonalizer for commutants of given finite group representations. As an application, when the given representation is orthogonal, we compute ordinary $QR$ factorizations blockwise in the isotypic basis. After transforming back to the original coordinates this yields a factorization $A=QX$, where $Q$ is orthogonal and $X$ is generally not upper triangular. Moreover, given a Wedderburn-adapted orthogonal basis, the transformed-back factors can be computed within the commutant.

We consider an optimal control problem governed by a stationary linear transport equation in two or three spatial dimensions, motivated by photon transport in radiotherapy treatment planning. The control is an inflow boundary source, and the observation operator may track either the angular flux or its angular integral, the latter serving as a simplified dose surrogate. The associated optimality system consists of a transport–adjoint pair with angular dependence, coupled through a boundary optimality condition. We formulate a residual least-squares functional whose minimisers coincide with solutions of the Karush–Kuhn–Tucker system and show that the functional controls the $L^2$ errors in the state and adjoint and the flux-weighted $L^2$ error in the control. This yields a continuum-level a posteriori estimate without recourse to discretisation. An additional unweighted optimality residual certifies the control in the norm of the objective. The estimate holds for every admissible approximation and accounts for the different boundary measures in the control penalty and transport traces, without a positive lower bound on $|\omega\cdot n|$. We then parameterise the state, adjoint and control using neural networks and minimise over a neural network space using appropriate quadrature. The Monte Carlo formulation replaces each residual integral by a weighted sample mean, whose quadrature error enters the a posteriori estimate. We prove convergence under explicit assumptions on approximation, quadrature and optimisation and derive bounds on the total error in terms of approximation, quadrature and optimisation error. Numerical examples in two and three dimensions illustrate the residual–error relation and the flexibility of the residual framework for field- and dose-tracking control problems, including regular manufactured solutions.

Fractional Brownian motion provides a useful framework for modeling random phenomena with anti-persistent temporal correlations that arise in applications such as anomalous diffusion, hydrology, finance, and geophysical processes. In this paper, we study a class of semilinear stochastic partial differential equations driven by additive fractional Brownian motion with Hurst parameter $H\in(0,\frac12)$. We establish the well-posedness and space-time regularity of the mild solution and investigate the strong convergence of a stochastic exponential integrator for the temporal discretization. The analysis allows the linear operator to be non-self-adjoint and combines analytic semigroup estimates with the canonical Hilbert-space structure associated with fractional Brownian motion and Malliavin calculus to handle the low temporal regularity of the noise. We derive an $H$-dependent strong temporal convergence rate, which reaches $H+\frac12$ under maximal spatial regularity. Numerical experiments for a stochastic advection–diffusion–reaction problem with heterogeneous Darcy flow confirm the theoretical temporal convergence rates. We also estimate the mean of the solution with different Hurst parameters $H\in(0,\frac12)$ using the Monte Carlo method.

Streaming GPU solvers for entropic optimal transport (EOT), such as FlashSinkhorn, avoid storing the dense kernel but still evaluate all $n\times m$ point pairs in every Sinkhorn iteration. We present FlashSinkhorn~2 (FS2), a solver for squared-Euclidean cost on low-dimensional point clouds that solves large discrete EOT problems to a prescribed marginal residual on a single GPU by coupling two stages. A coarse stage solves on cell centroids, lifts the potentials to every point and, when a sampled marginal check rejects the lift, continues on the centroids, replacing most point-level updates. A block-sparse fine stage then removes the centroid error that coarse updates cannot. Its Morton-ordered blocks support screening and fused tensor-core execution, and a threshold set by the block masses bounds each omitted tile's contribution to every row and column. On synthetic benchmarks, FS2 reaches the target residual on all 32 problems and GeomLoss multiscale on 10. On one A100, FS2 solves discrete EOT between two $1.34\times10^8$-particle measures from a cosmological $N$-body simulation, at an entropic blur equal to the mean interparticle distance, to an all-particle marginal residual below 0.01 in under 2.5 hours. To our knowledge, it is the largest discrete EOT problem solved to this accuracy within hours. For reproducibility, we release an open-source implementation at https://github.com/ot-triton-lab/flash-sinkhorn

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 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 of the underlying approximate posterior in total-variation distance and Kullback-Leibler divergence and, moreover, propose a hybrid sampler that is asymptotically exact by Metropolization. We validate the accuracy estimates and demonstrate the amortization in several numerical examples, including electrical impedance tomography and a likelihood-free Lotka-Volterra model.

We develop a finite element method (FEM) for partial differential equation (PDE) solver based on the entanglement-basis representation introduced in our companion work. By lifting non-linear finite-element equations into an augmented coefficient space, the governing PDE together with boundary, initial, and inter-element constraints can be expressed through a unified quadratic residual minimization. Although this augmented space grows exponentially with the number of elements, its tensor-product structure allows it to be represented efficiently using tensor networks. Using the matrix product state (MPS) as a concrete example, we show that density matrix renormalization group (DMRG) sweeps enable element-by-element optimization without explicitly constructing the full augmented space. For bounded bond dimension, the resulting computational cost scales polynomially with the number of finite elements. We extend the framework to time-dependent problems through implicit temporal discretization and demonstrate convergence under both mesh and polynomial refinement using diffusion equations.

We study free-boundary problems within a physics-informed framework using Kolmogorov-Arnold network (KAN) approximations. The proposed approach incorporates obstacle constraints, partial differential equation (PDE) inequalities, complementarity conditions, and boundary conditions through residual-based loss functions. We consider a linear elliptic obstacle problem, a nonlinear $p$-Laplacian obstacle problem, and a time-dependent one-phase Stefan problem. The proposed KAN solver is compared with physics-informed neural network (PINN) and residual-network baselines. Numerical experiments show that KANs achieve low relative $L^2$ and $L^\infty$ errors while accurately resolving contact regions and moving interfaces. The results indicate that KAN representations provide an effective alternative for solving free-boundary PDEs.