scientific computing teams

Abstract: Tensor decoupling algorithms recover compact representations of multivariate functions by combining linear factor matrices with univariate internal nonlinearities. We identify two limitations of recent methods. First, their performance is not invariant to mismatched output or gradient magnitudes, and second, their spline internals require manual tuning of their degrees of freedom. We propose two lightweight fixes. The first is a preprocessing step that rescales the input Jacobian tensor, removing both the bias towards larger magnitude outputs and the ill-conditioning of certain factor updates caused by such a mismatch. The second replaces B-spline internals by penalized splines coupled with generalized cross-validation, making the method less dependent on the user selecting the right number of control points. Finally, we release an efficient open-source JAX implementation of our algorithm.

Abstract: Implicit time integration is the standard way to overcome the tight CFL time step restriction when solving nonstationary PDEs numerically. However, within an adaptive time accurate integration, the time step size selection is governed not only by stability but also by accuracy requirements. Therefore, an eligible question to ask is: is (unconditionally) stable implicit time stepping really needed if the time step is restricted by accuracy considerations anyway? And how much do we then lose (in terms of the time step size) when switching from the expensive implicit time stepping to stabilized explicit schemes? An answer to this question is also important because such a switch is often done when solving large scale nonstationary PDEs on massively parallel computers. This paper attempts to answer these questions by numerical tests for two particular adaptive time accurate settings: (i)~nonstationary convection-diffusion equations and (ii)~nonlinear heat equations.

Abstract: Numerical libraries can still lose many bits of accuracy on localized ill-conditioned inputs, such as roots, poles, cancellation points, and branch boundaries. This paper presents NPTaylor, a Taylor-approximation-based framework for repairing localized high floating-point errors in numerical C routines. Given a high-error input, a localized repair domain, a high-precision oracle, and an error threshold, NPTaylor searches over Taylor degree, expansion center, segmentation, and floating-point evaluation strategy. It emits guarded C patches that use a generated Taylor evaluator inside the repair domain and fall back to the original routine outside it. Candidate patches are accepted only after validating the emitted C code on recorded local samples, guard endpoints, and the pre-repair peak input. The generated evaluator uses coefficient splitting and compensated Horner evaluation to reduce rounding error introduced by the patch itself. We evaluate NPTaylor on GSL benchmark workloads. On 38 one-input functions, both NPTaylor and AutoRNP repair all 114 configurations across three accuracy thresholds. NPTaylor is 3.33x faster by geometric-mean offline repair time and emits 0.162x as many C bytes. In the strictest setting, its median repair time is 3.74s versus 81.97s for AutoRNP. NPTaylor also repairs localized multi-input benchmark cases beyond AutoRNP's implemented one-input workflow. These results show that Taylor-based guarded patches can provide compact, efficient local repairs, while the claims remain validation-relative rather than full-domain proofs.

Abstract: Many machine learning approaches are two-level: an expert selects a class of methods, within this class the choice of key components is critical and usually manual. This paper addresses the automatic tuning of these components and proposes a framework transferable across models of different nature. The multi-agent large language model (LLM) framework has three roles---analyst, coder, and recorder~---with an experiment log and a sandboxed execution environment. Unlike systems that optimize the entire ML script or search architectures in a discrete space, the framework operates over an open but semantically narrow space: the agent generates code only for the specified components, leaving the overall network structure unchanged. The approach is validated in two directions: tuning a classification spiking neural network with correlation learning, and solving high-order differential equations with Physics-Informed Neural Networks (PINN). Across five tasks---spiking classification on MNIST and Digits, and boundary-value problems for the generalized nonlinear Schrödinger equation (GNLSE) and Kuramoto--Sivashinsky (KS) equations---the agent outperforms the baselines: classification accuracy is 13 percentage points (MNIST) and 6 percentage points (Digits) higher; the error against the analytical single-soliton solution is reduced approximately 2491-fold; the invariant-preservation error drops approximately 325-fold in the selected GNLSE run with a Gaussian perturbation and by more than 100-fold on average over repeated initializations; and the mass drift drops 14-fold for the KS problem with a Gaussian perturbation. The framework thus tunes key components on tasks of different nature without changing its internal structure.

Abstract: Large-scale scientific simulations generate substantial data volumes, making lossy compression essential for reducing storage and data movement costs. However, users configure compressors through numerical error bounds (EBs) while often evaluating results using quality metrics and end-to-end performance. Because the relationship between an EB and these outcomes varies across datasets and compressors, identifying a suitable configuration typically requires exhaustive search, which can be time-consuming and computationally demanding. We present LOCAA, an Large language model-based scientific lossy compression auto-tuning agent that performs compression-in-the-loop search using tool-integrated execution, compressor-aware guidance, and persistent memory. LOCAA supports user-defined objectives and constraints without requiring a specialized search strategy. We evaluate LOCAA across three scientific applications, five compressors, and three use cases: fixed-ratio compression, compression tuning under multiple quality constraints, and compression tuning under quality and time constraints. For fixed-ratio search across twelve fields from two applications, LOCAA requires 1.98x fewer evaluation trials than binary search and 5.03x fewer than FRaZ on average. For compression ratio maximization under joint Peak Signal-to-Noise Ratio (PSNR) and Structural Similarity Index Measure constraints across six fields from the Community Earth System Model application, LOCAA reduces the average evaluation trials from 61.5 using binary search to 17, corresponding to a 72.4% reduction. Persistent memory further reduces the average number of trials by 27.9% across multiple timesteps of the same field. These results demonstrate the potential of tool-augmented LLM agents to provide flexible and efficient compressor tuning across diverse scientific data and user-defined objectives.

Abstract: Agentic artificial intelligence systems can carry out a broad range of tasks in software engineering and scientific research, from writing and translating code to running workflows for data analysis and visualization. In scientific computing, the difficulty is verifying that agent-generated code is both correct and understandable to teams whose members bring different areas of expertise. We therefore argue that these systems are best used within collaborative team structures rather than as full automation. In the workflows we propose, domain experts write the specification and plan, and agents operate under a deterministic orchestration pattern to write the target code. Each stage ends with a numerical comparison against the reference code and requires human review and approval before the next begins. We evaluate these workflows on the translation of a large high-energy physics application from Fortran to C++, running the same task under different orchestrators, design patterns, and models. Across fourteen experiments, a simple author--reviewer loop with enforced limits completed a comparable number of files to a multi-agent workflow at about one-third of the cost per file.

Abstract: We present three practical tutorials on numerical computation and machine learning for mathematical research, developed for the DANGER: Data, Numbers, and Geometry workshop held at the Banff International Research Station in April 2026. The first develops a numerical approach to exterior calculus from pointwise evaluations of differential forms, using a flux formulation of the exterior derivative. Examples in Euclidean space and on the sphere illustrate geometric identities, topological features, and the effects of approximation and finite precision. The second examines how mathematical structure guides neural network design through examples involving elliptic curves, quivers, and a boundary value problem. It explores how architectural choices affect learning and uses interval arithmetic to bound the residual of a trained network over the full interval of the boundary value problem. The third addresses the evaluation and presentation of machine learning results, covering performance metrics, statistical uncertainty, classification thresholds, receiver operating characteristic curves, and accessible figure design. Throughout, the tutorials distinguish numerical agreement, predictive accuracy, structural guarantees, and rigorous bounds as different forms of evidence. Each contribution can be read independently, with accompanying notebooks and exercises that allow readers to reproduce the examples and adapt the methods to other problems.

Abstract: Ensemble guidance combines pretrained diffusion priors with black-box forward models to solve inverse problems without differentiating through the physical simulator. However, observation coordinates with large predictive spread or extreme residuals can dominate the ensemble correction, degrading reconstruction accuracy. We show that two simple modifications, weighting and clipping, substantially improve this correction. Our method, Robust Ensemble Guidance (REG), uses ensemble predictive spread to balance observation scales and adaptively clips standardized residuals to limit the influence of extreme discrepancies. Both operations reuse existing particles and forward predictions, requiring no additional denoiser or forward-model evaluations. Under a local linear Gaussian model, we derive conditions for reduced one-step estimation risk, bound the influence of individual observation coordinates, and characterize when these benefits persist with finite ensembles. Experiments on Navier-Stokes inversion, black-hole imaging, and acoustic full-waveform inversion demonstrate improved reconstruction over the underlying ensemble solver. In particular, REG increases black-hole reconstruction PSNR by 6.2-8.2 dB across three observation regimes and reduces Navier-Stokes reconstruction error by 26.4\% in a matched-budget comparison. These findings highlight the importance of observation heterogeneity and residual influence in designing reliable generative solvers for scientific inverse problems.

Abstract: Current and next-generation large-scale experiments, especially in high-energy physics and cosmology, involve increasing volumes of data that require exascale computing power to analyze and model. Exascale machines present challenges for workflow scaling and efficient resource usage. We propose a multi-stage task scheduling approach for optimizing resource usage for these workflows. We demonstrate task scheduling for an example Monte Carlo simulation pipeline from the Short-Baseline Near Detector (SBND), and present theoretical and numerical tools for modeling task scheduling in these workflows that can be used to identify optimal scheduling strategies and choose resources for a given problem size.

Abstract: Implicit neural representations (INRs) provide a flexible framework for modeling high-dimensional continuous fields, but their training is often inefficient due to uniform subsampling that ignores spatial heterogeneity. Existing adaptive sampling methods partially address this issue by prioritizing high-error samples, but typically operate at the point level, often leading to redundant sampling in localized regions and insufficient coverage of the domain. We propose ACES (Adaptive Coverage-aware Efficient Sampling), a structured sampling framework that improves training efficiency by decoupling coverage and importance. ACES constructs adaptive spatial partitions to ensure domain coverage and reduce redundancy, and applies region-level importance weighting to prioritize informative regions during training. We provide a theoretical analysis showing that adaptive partitioning reduces gradient variance by increasing within-region homogeneity, and that controlled bias in region-level weighting may improve optimization efficiency relative to standard unbiased estimators. Experiments on scientific field learning tasks demonstrate that ACES achieves faster convergence and lower error than uniform and pointwise adaptive sampling baselines, with the largest gains in fields with highly localized complexity.

Abstract: Diffusion models are a powerful generative paradigm used across multimedia and scientific applications. Guided diffusion methods impose requirements on the generation by adding the gradient of a differentiable loss (the guidance function) as a drift term during inference. The weight of this drift (the guidance scale) is critical for the trade-off between data quality and requirement satisfaction. To achieve both of these goals, guided diffusion must resort to small guidance scales and lengthy sampling, incurring high computational costs. This work proposes learned end-to-end guidance schedules (LEEGS) to achieve these objectives with fewer sampling steps. LEEGS trains a time-dependent schedule by minimizing the guidance function over a small set of examples using stochastic gradient descent. Backpropagating through guided sampling is computationally expensive, so LEEGS uses an approximation of the gradient that cuts training time by a factor of 4. We evaluate LEEGS on diverse guidance tasks, including (a) image inpainting, (b) noisy image inverse problems, (c) face-ID-guided generation, and (d) forward and inverse PDE problems, outperforming baselines at equal budget (50 or 100 NFEs), or matching constant guidance with only 10% of the steps.

Abstract: Many forms of data, including physical fields, geometric shapes, and visual signals, are naturally described by functions over continuous domains but are observed through discrete samples. Representing these functions on fixed uniform grids imposes a trade-off between resolving localized variation and increasing computation across the domain. Neural operators address this mismatch by learning mappings between functions, while latent-attention architectures provide flexible processing of sampled observations. We introduce the Function-Space Transformer (FST), a framework for learning from functions through a spatially adaptive continuous latent representation. FST stores features at anchors whose locations are predicted from the input observations and recursively refines these anchor features through function-space interactions. This allows the representation to adapt its spatial organization to each input rather than inherit that of the observation grid, while supporting both spatially resolved and finite-dimensional outputs. On PDE solution prediction using PDEBench Burgers and Darcy flow, FST substantially outperforms the Perceiver IO baseline, whose latent representation lacks explicit spatial organization, and is highly competitive with the Fourier Neural Operator. On ImageNet-1K, FST achieves higher classification accuracy than the Vision Transformer baseline, with fewer parameters across these comparisons. Ablations further support the benefits of function-space updates and recursive refinement. Together, these results highlight the potential of adaptive continuous representations for both scientific prediction and visual recognition.

Abstract: Reduced and mixed precision have moved from a niche optimization to a central design axis in scientific computing and engineering, driven by energy constraints, heterogeneous accelerators, and the convergence of simulation and machine learning. This paper organizes the landscape around seven coupled themes---number formats, floating-point emulation, emerging architectures, hardware/software co-design, relation to other approximations, software design, and precision as a multilevel resource ---and, for each theme, synthesizes the state of the art, future directions, and open questions. We emphasize \emph{energy per trusted solution} as the core objective, and we frame \say{recklessly responsible} computing as a pragmatic doctrine: exploit low precision aggressively, but with systematic detection, escalation, and certification pathways.

Abstract: The Girard--Hutchinson estimator provides an extremely simple randomized estimate of the Frobenius norm of a matrix $A$ that can only be accessed implicitly via matrix-vector products. In particular, if $Ω$ is a random Gaussian matrix with $r = O(1/\varepsilon^2)$ columns, than $\frac{1}{r}\|AΩ\|_F^2$ provides a $(1\pm \varepsilon)$ multiplicative approximation to $\|A\|_F^2$ with high probability. In this work, we introduce a closely related estimator, given by \begin{align*} {\frac{1}{r}\|AΩ\|_F^2 + \frac{1}{r}\|Ψ^T A\|_F^2 - \frac{1}{r^2}\|Ψ^T AΩ\|_F^2}, \end{align*} where $Ψ$ is a second, independent random Gaussian matrix with $r$ columns. We prove that this estimator yields a $(1\pm\varepsilon)$ multiplicative approximation to $\|A\|_F^2$ when $r = O(1/\varepsilon)$, a quadratic improvement over Girard--Hutchinson. This dependence on $\varepsilon$ is optimal. Our method, which we call Hutch# (pronounced ``Hutch sharp''), matches the complexity of the Hutch++ algorithm [Meyer, Musco, Musco, Woodruff, 2021]. However, unlike Hutch++, Hutch# uses only \textit{non-adaptive} matrix-vector products with $A$ and $A^T$ and requires no orthogonalization or other adaptive linear algebra steps. Thus, Hutch# combines the simplicity of the Girard--Hutchinson estimator and the optimal query complexity of Hutch++.

Abstract: We evaluate a SymPy implementation of the parallel Risch-Norman method using Charlwood's 2008 suite of 50 challenging indefinite integrals. The system automatically builds integrand towers, verifies answers by differentiation, and returns correct, verified integrals for 49 problems with zero errors. The single failure, $\int\arcsin(x\sqrt{1-x^2})\,dx$, is proven non-elementary using a holomorphic-remainder certificate, though limited by an unverified completeness hypothesis on a genus-three curve. Part II's degree bounds successfully reduce classical ansatz sizes by two-thirds without impacting running time. The paper details algorithmic mechanisms like $S'$-units, Pell units, and residue computing in tower coordinates. Compared to mature implementations, this untuned SymPy prototype is slower than FriCAS (by a factor of six on the median integral) but more accurate than AXIOM (which returned three wrong answers). Profiling pinpoints performance bottlenecks in nonlinear norm searches and nested number fields, outlining clear targets for optimisation.

Abstract: Physics-informed neural networks (PINNs) use one shared representation over the computational domain, which can become difficult to optimize on long domains and for high-order operators. We study a minimal alternative: multiply the first hidden activation of an otherwise unchanged dense PINN by input-dependent localization functions, giving first-layer units receptive fields without partitioning the domain or adding interface losses. We screen 13 families of localization functions, in up to three parameterizations each, on a nonlinear harmonic oscillator (HO), a heat equation on a long spatial interval, and a manufactured four-dimensional (4D) fourth-order problem, with ten paired seeds throughout. Three configurations give large reductions in solution error at matched budgets: (i) Fixed Gaussian localization functions on the $2π$ HO domain cut mean solution RMSE from $4.8369\times10^{-1}$ to $8.83\times10^{-3}$ at 3k epochs. (ii) The inverse-quadratic family with learnable centers and widths cuts it from $2.896\times10^{-1}$ to $3.06\times10^{-2}$ on the $8π$ heat domain at 10k epochs. (iii) Fixed bump localization functions cut it from $1.75947\times10^{1}$ to $2.260\times10^{-1}$ on the $4π$ 4D domain at 10k epochs. Every paired seed improves in these three comparisons. The screen also shows that the mechanism is not a free win: on HO only 2 of 13 families beat the baseline, and 10 of the remaining 11 are 9 to 23 times worse; on 4D four families are non-finite and five are more than three orders of magnitude worse than the baseline. The inverse-quadratic family is the only one that beats the baseline on all three equations. Overall, these results show that first-layer localization can provide measurable improvements to baseline PINNs on long-domain and high-order problems.

Abstract: Sparse computations are important workloads in applications such as scientific computing, graph neural networks (GNNs), and machine learning. While many sparse operations can benefit from modern GPUs, the sparsity pattern remains important to performance because it affects memory coalescing, block organization, and load balancing. Previous studies show that spectral signatures can help analyze the global structure of sparse matrices. The fast Fourier transform (FFT) is commonly used to extract spectral signatures, and efficient GPU FFT libraries are available. However, sparse matrices, especially adjacency matrices for large graphs, tend to be very large and sparse. Existing dense-matrix-based FFT implementations are difficult to scale up, making the spectral patterns of these matrices difficult to obtain. We therefore propose a three-fold research approach comprising a lossless Binary-Sparse FFT (BS-FFT) and two compression methods: Elastic BS-FFT, which reuses the BS-FFT pipeline on a sampled frequency grid, and density-map-based spatial compression. Experiments show that BS-FFT reduces GPU memory use by 2.9--11.6 times relative to dense cuFFT and completes all 15 GNN adjacency matrices where dense cuFFT completes 6 on a 40 GB A100. Elastic BS-FFT and Density Map compression reduce GPU computation time by 2.0--1466.4 times relative to BS-FFT with spectral feature error of only 0.16% to 11.56% across the sampling rates from 6.25% to 0.0061%.

Abstract: We consider the problem of estimating the trace of an implicit matrix $\mathbf{A} \in \mathbb{R}^{d^p\times d^p}$ that can only be accessed through matrix-vector products queries. The \textit{Hutchinson trace estimator}% ~\cite{Girard1987algorithme, article-hutchinson} is a classical sketching method for this problem. Their estimator, $H_{m}(\mathbf{A}) = \frac{1}{m} \sum_{i=1}^{m} {\mathbf{z}^{(i)}}^T \mathbf{A} \mathbf{z}^{(i)}, \quad \text{where } \ {\mathbf{z}^{(i)}}\in \mathbb{R}^{d^p}$, and $z^{(i)}_j \in {N}(0, 1), j\in [d^p]$, satisfies the following guarantees: (i) $\mathbb{E}[H_{m}(\mathbf{A})]=\operatorname{tr}(\mathbf{A})$, and (ii) $\mathrm{Var}[H_{m}(\mathbf{A})]=\frac{2}{m}||\mathbf{A}||_F^2$. Generating one query vector $\mathbf{z}^{(i)}$ requires $O(d^p)$ random bits; thus, $m$ queries require $O(md^p)$ random bits, which can be prohibitive in large-scale applications. Recent work by Meyer et al.~\cite{meyer2025hutchinsonsestimatorbadkroneckertraceestimation} proposes a variant of the Hutchinson trace estimator in which each query vector in $\mathbb{R}^{d^p}$ is constructed as the Kronecker product of $p$ random vectors in $\mathbb{R}^d$, requiring $O(mpd)$ random bits for $m$ query vectors. The estimator of~\cite{meyer2025hutchinsonsestimatorbadkroneckertraceestimation} is unbiased; however, its variance grows exponentially with $p$. In this work, we address this limitation by proposing a sketching-based estimator that requires $O\!\big(p (d + m)\log m\big)$ random bits, yields an unbiased estimate of the trace, and simultaneously achieves a variance bound that grows polynomially with $p$.

Abstract: Scientific applications increasingly rely on high-performance computing (HPC), yet translating a scientist's high-level goal into a correct target-scale execution remains brittle and labor-intensive. Large language model (LLM) agents promise to automate this, but two obstacles remain: granting a cloud-hosted model direct HPC access exposes credentials and execution authority, while withholding it demands continuous human supervision; and one-shot generation cannot adapt when generated artifacts fail in a site-specific HPC environment. We present \textsc{ECAS}, an \textbf{E}dge-\textbf{C}ontrolled \textbf{A}gentic \textbf{S}ystem for closed-loop execution of scientific computing campaigns with limited human intervention. \textsc{ECAS} separates \emph{reasoning}, \emph{control}, and \emph{execution}: a cloud-hosted LLM proposes plans, artifacts, and repairs; a user-controlled edge agent retains credentials, workflow state, and execution authority while enforcing policy and resource constraints; and the HPC system computes. Its core mechanism is \emph{validation-gated execution}: generated artifacts pass static checks and small-scale validation, failures trigger repairs from sanitized execution feedback, and target-scale execution is permitted only after validation and policy checks pass. \textsc{ECAS} also draws on an edge-resident library of expert-distilled, site-specific skills that is never disclosed to the cloud. In preliminary experiments with three scientific applications on two production ALCF systems under six injected fault types, closed-loop repair improves application success from 0/6 to 6/6 over one-shot generation, validation gating prevents all three observed target-scale failures, and skill conditioning improves success from 4/6 to 6/6. These results show the feasibility of delegating adaptive reasoning to the cloud while retaining execution control at the edge.