numerical software developers

Abstract: Solutions of second-order elliptic PDE problems on nonconvex or cracked domains with mixed Dirichlet-Neumann conditions typically lie outside $H^2(Ω)$. Near re-entrant corners, crack tips and points where the boundary condition changes type, the gradient grows like $r^{λ-1}$ with $λ<1$. We study residual-based adaptive refinement in hybrid variational-collocation PINNs. A multiscale Fourier-feature network is trained on a composite loss. The loss combines a variational residual tested against Legendre bubble functions, a strong-form collocation term and penalties for the Dirichlet and Neumann data. The variational term needs only first derivatives of the network, whereas the collocation term requires second derivatives and provides a residual that can be evaluated at any point. After each training stage, candidate points selected by this residual are added to the collocation set. The method is tested on five two dimensional benchmarks, namely a Mode~III crack, a cross-shaped domain with four re-entrant corners with unequal amplitude, a reaction-diffusion problem with four localized peaks, a monotone semilinear problem with a boundary layer and the Motz problem. Compared with uniform point insertion under the same collocation and iteration budget, residual-based refinement reduced the relative $L^2$ error by up to a factor of $32$, the relative $L^\infty$ error by up to a factor of $38$ and the relative $H^1$ semi-norm error by up to a factor of $15$.

Abstract: We study how many bits are needed to write down an elementary analytic expression that uniformly approximates a Lipschitz function. Expressions are finite trees over the basis $\{+,-,\times,÷,\exp\}$ with the variables and rational constants, and their size is $\dc(f)$, the length of the shortest description of $f$ in a fixed binary description language. For functions of one variable we prove that every $L$-Lipschitz $F:[0,1]\to[0,1]$ has an everywhere-defined approximant of accuracy $Δ$ with $\dc(f)\le C\,(L/Δ+\lg(1/Δ+2))$ for all $L,Δ>0$, by a linear spline smoothed with logistic ramps in a description language with one repetition construct. We prove a matching lower bound valid for \emph{every} description language: for $Δ\le 1/36$ and every $L>0$ some $L$-Lipschitz function needs at least $\tfrac1{16}(L/Δ+\lg(1/Δ+2))$ bits. Hence the worst-case description complexity is $Θ(L/Δ+\lg(1/Δ+2))$. For functions on $[0,1]^d$ the same method gives $Θ_d\bigl((L/Δ)^d+\lg(1/Δ+2)\bigr)$, using a tensor-product spline whose coefficients are mixed second differences and cost $O(d)$ bits each. Without the repetition construct the upper bounds carry an additional logarithmic factor, whose necessity for plain expression trees is left open. A reference encoder, decoder and evaluator confirm the constructions numerically, and figures and tables report the checks.

Abstract: Neural PDE surrogates are trained on numerical solver outputs that contain both physical evolution and solver-specific discretization errors. Because surrogates are also evaluated against held-out trajectories from the same solver, standard benchmarks cannot distinguish fidelity to the exact evolution from imitation of the numerical scheme. We introduce an empirical Fourier-symbol diagnostic that probes a trained surrogate's linearized one-step operator with individual Fourier modes and compares it with both exact-evolution and training-scheme references. To address architectural spectral bias, we train identical networks on schemes with orthogonal dissipative and dispersive signatures and compare their learned operators. In linear advection, the learned surrogates reproduce the training schemes' amplitude and phase errors, with the twin-scheme difference reaching more than 99.8\% of the analytically predicted full-imitation ceiling. The same behavior occurs for a non-local Fourier neural operator and at the operator level for nonlinear Burgers dynamics. These results show that agreement with solver-generated test data does not by itself establish fidelity to the exact evolution. Fourier-symbol measurements provide a direct diagnostic of numerical provenance.

Abstract: We prove that the probability that Gaussian elimination with partial pivoting on $n \times n$ random matrices has growth $ρ$ is at least inverse quasi-polynomial: $Ω(\exp(-c \log^2 (ρ)\log(n)))$ for some constant $c > 0$. This lower bound breaks standard conjectures in the smoothed and average-case analysis of Gaussian elimination. To the best of our knowledge, it is the first non-trivial lower bound on the probability of large growth for Gaussian random matrices.

Abstract: 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.

Abstract: We present a novel algorithm for calculating the preimage of an affine space through a product ${J}={J}_{T-1}\cdots{J}_0$ of matrices ${J}_t$ of special form: finding the largest input space $\mathbf{X}$ such that $\mathbf{x}\in\mathbf{X}$ implies ${J}\mathbf{x}\in \mathbf{Y}$, where $\mathbf{Y}$ is a given output affine space. These special matrices arise in AD, where the Jacobians $J$ describing the linearized computation have precisely this structure: the product of a series of linearized primitive numeric operations. This allows us to use the new algorithm to formulate Null-A mode preimage AD, which finds the affine preimage through the Jacobian or Jacobian transpose of a numeric computation. This is a generalization of the inverse AD problem of solving ${J}\acute{\mathbf{x}}^{\ast}=\acute{\mathbf{y}}^{\ast}$ or ${J}^{T}\grave{\mathbf{y}}^{\ast}=\grave{\mathbf{x}}^{\ast}$. The key is to represent affine spaces in a fashion which lends itself to efficient preimage calculation, in a compositional and quasi-local fashion, through a succession of matrices ${J}_t$. Unlike previous methods, Null-A preimage mode AD allows the ${J}_t$ matrices to be non-square, corresponding to a computer program whose number of active variables swells and shrinks during the computation. When ${J}$ is square and the initial affine space is a single point, this finds the conventional inverse. But in the more general case, having the entire affine space provides freedom which can be leveraged in a problem-specific manner. We apply the method to small problems on-CPU where the $J_t$ are linearized scalar unary or binary numeric functions; and to larger problems on-GPU where the $J_t$ are linearized aggregate array operations like convolution and attention.