Recursive Improvement of a Differentiable Scientific Software Ecosystem
Authors: Pengcheng Hou, Xiaojun Tan, Sihan Hu, Ruisi Wang, Shuo Chen, Lei Wang, Youjin Deng, Kun Chen
Organizations: Hefei National Laboratory, Hefei 230088, China · Institute for Advanced Algorithms Research, Shanghai 200120, China · University of Science and Technology of China, Hefei, Anhui 230026, China · Institute of Physics, Chinese Academy of Sciences, Beijing 100190, China · Institute of Theoretical Physics, Chinese Academy of Sciences, Beijing 100190, China
Differentiable programming connects scientific computation with gradient-based inference, learning and design. Extending these capabilities across a heterogeneous software ecosystem requires specialized effort to implement derivatives, integrate interfaces and evaluate quality. AI coding agents can accelerate this transformation, but translating their capabilities into useful scientific software requires identifying research needs and evaluating how well implementations meet them. We present an environment for agent-driven evolution of differentiable scientific software that connects demand identification, development and quality evaluation. A unified differentiation interface exposes reusable derivative rules alongside existing numerical routines, allowing research tasks to share these capabilities. Research requirements guide development, with implementations assessed through independent derivative checks, workflow tests and performance evaluation. Validated software, research programs and tests become shared resources for subsequent studies. We construct and validate automatic differentiation extensions across 20 packages spanning physical, chemical and biological modeling, with research workflows demonstrating reuse across tasks. Benchmarks demonstrate computational savings over finite differences in gradient evaluation and complete parameter estimation. Research-driven revisions make previously unsupported workflows differentiable, correct derivatives of scientific observables and eliminate redundant computation. Quantum-control and thermal-design studies revise objectives in response to physical evaluation, improving designs while reusing existing derivatives. This work provides a practical approach to expanding differentiable programming across established scientific software and organizing AI agents around the recursive improvement of a shared computational ecosystem.
Figures & tables
Figure 1: AD extensions and recursive improvement of a scientific software ecosystem. (a) Extensions connect established computations to gradient-based research across 20 packages, using runtime-specific value-and-derivative interfaces. (b) Research reuses shared resources and returns programs, checks and development requirements. Updated software and tests then support further research. Papers, public research code and LKM inform task discovery (Section 2.2 ).
Numerical capability
Research computations using it
Validation evidence
Diffusion sensitivities FiPy
Temperature-based conductivity inversion; thermal-cloak optimization under different heating objectives.
Native field replay and finite-difference gradient checks in both studies (Appendices E , F ).
State-propagation derivatives QuSpin
Nominal and coupling-aware state-preparation optimizations: one propagation rule, 96 controls.
Independent matrix-exponential replay and directional finite differences (Appendix D ).
Evolution, steady-state and spectral rules QuantumOptics
Value and gradient checks of composed control, inference and spectral objectives.
Table 1: Research uses and checks of reusable derivative capabilities. Rows distinguish optimization studies, composed-workflow checks and sensitivity cases. Workflow counts are reported evaluations; definitions and sources appear in Appendices B – C .
Figure 2: Computational utility of derivative implementations across 20 scientific packages. Each row measures a specified gradient workload. The component column gives P/D : FD evaluates P derivatives, while AD returns a gradient over all D independent parameters. Solid bars have P=D ; hatched bars have P<D . Bars show ratios of median runtimes, tFD(P)/tAD(D) , with absolute times below. † : reimplemented diffusion kernels checked against original-software outputs. MCell compares an analytic expectation rule with native Monte Carlo FD under a 10% gradient-accuracy criterion. Workload definitions and repeats appear in Table 6 ; validation details appear in Appendices B.4 and B.6 .
Figure 3: A research task drives a new derivative capability and then an implementation improvement in Kwant. (a) Eight-gate potential designed using the first scattering-rule update. (b) Native transmission of initial and AD-designed devices at 101 energies; circles mark four reserved points per energy window. Dashed thresholds apply to their mean transmission within each window. Two incoming channels give a maximum of two. (c) Median time and all three process measurements for one value and full-gradient request at the same saved design, with matching numerical results; the gray line is the median FD reference for the same request. The feedback band connects the missing derivative, its use in design and the later execution improvement.
Figure 4: A revised control objective improves robustness using the same QuSpin derivatives; saved states then support a measurement study. (a) Gradients pass through evolution to 96 controls. The pulse map shows the selected robust control. Readout probabilities use coupling shift δ=0.1 and encoded phase ϕ=π/32 . (b) Error 1−F for three paired initializations, where F is GHZ squared overlap. Line styles identify pairs; squares mark optimization conditions and circles other evaluations. (c) A subsequent sensing task reuses two selected, unpaired states. Phase-estimation root-mean-square error (RMSE) is evaluated at known δ=0.1 and ideal reference probes. Each point uses 500 repetitions of 10,000 measurements. Markers distinguish phases; bars show approximate Monte Carlo standard errors. Reference probes have no matched preparation cost. Protocols appear in Appendices D and D.2 .
Figure 5: A thermal cloak designed with FiPy diffusion derivatives following a joint objective and training-grid revision. The eight-layer design is optimized on a 192×192 grid and evaluated at 576×576 with material parameters fixed. (a) AD-optimized conductivities. (b) Cubic-load temperature deviations from the homogeneous reference, on a shared scale. Interiors are masked; dotted circles mark the observation annulus. (c) RMS disturbance relative to the bare core. The line separates three design loads from two new loads. Ratios below one (blue) indicate improvement; above one (orange), greater disturbance. Conductivity and ratio color scales are logarithmic. Both optimized designs use AD; the uniform coating is an analytic reference.
Appendix figures & tables17 assets
Supplementary material from the paper’s appendix.
Appendix
Stage
Inputs
Output used by the next stage
Investigation
Original software, generic development instructions, working runtime
Software-specific specification, numerical targets and runtime observations
Design review
Proposed specification and separate original-software checkout
Findings on scope, numerical semantics and differentiation route
Revision
Specification and review findings
Revised design and implementation targets
Implementation
Revised design, original code and numerical references
AD extension, documentation, package tests and explicit unresolved cases
Appendix
Table 2: Inputs and outputs of the four-stage development procedure. The reviewer receives the proposed design but inspects the original software independently. Each stage passes its documented outputs to the next.
Software
Evaluated path and original computation
Numerical check
FiPy
Spatial diffusivities and boundary inputs → native diffusion fields through implicit solve rules reused in inverse and cloak tasks
Native fields and FD gradients; boundary-map source tests
OpenPNM
Throat conductances → steady Stokes-flow pressure on a fixed pore network
Recorded solver-output and selected-coordinate FD checks
PHREEQC ( Parkhurst & Appelo, 2013 )
Solution inputs → weighted ionic strength, molalities and log activities
Real-engine speciation and all-input FD checks
QuantumOptics
Hermitian controls → observable of the original steady-state solver through a shared adjoint
Native steady solver and full-gradient FD check
QuSpin
Static Hamiltonian parameters → evolved state, control loss and observables through shared propagation rules
Independent propagation and FD checks; control evaluations
BioDynaMo
Voxel concentrations → one diffusion update through a reimplemented field kernel
Kernel/native-output parity and selected-coordinate FD
Appendix
Table 3: Evaluated parameter-to-output capabilities across 20 scientific software packages. Each row identifies the tested numerical boundary, its use of the original computation and its validation. These are evaluated paths, not complete API inventories. Citations identify the original software.
Software
Shared recorded interface
Cases
Paired
QuSpin
hamiltonian.evolve
10
2
OpenPNM
Transport run and flux rate (each)
66
66
TeNPy
dmrg_spectral_observable
9
9
TeNPy
tebd_spectral_observable
5
5
Phonopy
homogeneous_frequencies
4
4
LALSuite
SimInspiralChooseFDWaveform
16
16
Appendix
Table 4: Shared interfaces appearing in routine-case execution records. The last two columns count distinct case identifiers with any captured call and with a paired AD/FD timing for that interface, respectively. Repeated calls and runs count once. Counts do not certify the full scientific workflow or identify different underlying derivative rules.
Software
Scalar objective and active input
Headline FD coordinates
FiPy
Weighted steady diffusion–reaction solution; 808 face diffusivities on fixed 24×16 grid.
Seeded 256 faces.
OpenPNM
Weighted Stokes-flow pressure; 565 throat conductances on fixed 300-pore network.
Seeded 256 throats.
PHREEQC
Weighted ionic strength, molalities and log activities; six seawater totals plus pH, pe, alkalinity, temperature and pressure.
All 11 named inputs.
QuantumOptics
Retr(Oρss) with O=n/11+0.08(a+a†) and a controlled 12-dimensional Hamiltonian; 144 real Hermitian-basis controls.
All 144 active controls; base and active set vary across the smaller- P series.
QuSpin
Real target-state overlap after 32×32 symmetric-Hamiltonian evolution; 528 independent real matrix coefficients.
First 256 symmetric directions.
BioDynaMo
Weighted output of one periodic 163 diffusion step; 4,096 input concentrations in an extension kernel.
Seeded 256 voxels.
Appendix
Table 5: Summary of scalar objectives and active inputs for the 20 workloads in Figure 2 . Weights, observation operators and stated model conditions are fixed; the last column identifies the FD request. Timing and component counts appear in Table 6 .
Figure 6: AD/FD numerical checks for the workloads in Figure 7 . (a) FiPy relative L2 gradient differences. (b) Kwant maximum absolute component differences. Points show all three process records per size; lines connect their medians. Dashed lines mark the configured relative tolerance for FiPy and the absolute component of Kwant’s combined absolute/relative tolerance. Differences from FD are numerical comparisons, not errors against an exact derivative.
Software
Computation
P/D
n
FD (ms)
AD (ms)
FiPy
24 × 16 diffusion-reaction solution
256/808
7
1,670
10.8
OpenPNM
300-pore steady Stokes-flow pressure
256/565
9
419
2.22
PHREEQC
seawater equilibrium observables
11/11
5
58.9
3.21
QuantumOptics
12D steady-state observable
144/144
3
4,570
16.5
QuSpin
32 × 32 symmetric-Hamiltonian evolution
256/528
5
74.5
0.592
BioDynaMo
extension-kernel 16 3 diffusion step
256/4096
9
20
0.12
Appendix
Table 6: Definitions and reported timing medians for all workloads in Figure 2 . P/D gives derivatives requested by FD / the complete task-level gradient returned by AD (after contraction for MCell); when P<D , FD requests a selected subset contained in the full AD result, not the same complete gradient. n is the source-reported repeat count, not independent processes. Times are milliseconds and follow the Figure 2 ratio order. MCell ( P=D=1 ) uses five α=1 requests for timing medians; two α=0.8 requests validate accuracy and are excluded from those medians.
Figure 7: Complete-gradient cost as the parameter dimension grows. (a,b) Value plus full gradient: FiPy fixes its 32×32 mesh; Kwant grows to N=P+4 sites at energy 0.23. Open points are within-process medians of three timings; lines connect medians over three processes. (c) FD-normalized AD cost, c=(2P+1)tAD/tFD , using method-specific medians, without process pairing. Each measurement is one value-and-full-gradient request. Derivative checks are in Figure 6 .
Research requirement
Retained response
Subsequent use or checked result
Kwant: tune gates to filter electron energies
Scattering derivative, then shared linearization
Filter targets met; saved device reused to verify matching gradients at lower cost (Fig. 3 )
StochKit2: differentiate the requested population bins
Fixed-bin rules and reusable simulated trajectories
Numerical reconstruction agrees; later reuse preserves outputs with calls 4→2
OpenPNM: include reaction sources and directed transport
Source-term, directed-conductance and flux rules in the report-tested patch
Reaction–diffusion and advection–diffusion join the checked paths: 42→66 of 80 pass
BioDynaMo: preserve the original mechanics updates
Rule matching sequential updates and storage reorder
Eleven partial mechanics cases pass after repair
QuSpin: use shared propagation derivatives in research programs
Eight programs migrated; two more pass ( 23→25 /47); later control reuses the rule
QuantumOptics: compose further inference and control tasks
Nine programs using existing evolution, steady and spectral rules
Workflow retests: 35→44 of 49 pass; no new kernel implied
Appendix
Table 7: Research requirements, retained software changes and subsequent uses. QuSpin’s migrations and the OpenPNM, BioDynaMo and QuantumOptics outcomes are reported workflow retests; later QuSpin control directly uses the package rule. The surrounding sections give source identities and checks.
Figure 8: Native transport design after the scattering-rule addition. (a) The gate-induced potential reconstructed from the saved AD design and the eight coefficients for initial, AD and FD designs. The heatmap excludes the constant onsite term 4. (b) All eight reserved-energy transmissions; the dashed targets concern band means. Only measured energies are plotted, with no interpolation across the transition. Both designs meet the criterion; they are the last accepted iterates at the budgets in Table 8 .
Design
Stop mean
Pass mean
Accepted updates
Assemblies
Initial
1.99856
1.99998
—
—
Package VJP
0.26580
1.73687
40
3,078
Central FD
0.26751
1.74083
38
6,000
Appendix
Table 8: Reserved-energy performance of the fixed-device designs. Each band mean uses four points. AD stops at its update limit; FD stops at its assembly limit and retains the last accepted iterate.
Figure 9: FiPy execution comparison across ten research workloads. Open symbols show three separate process measurements per method and workload; filled symbols are medians. All ten workloads are retained, including Shape* with its unresolved derivative check and the four later-AD medians slower than FD. Revisions combine package preparation reuse and reverse-mode task code.
Figure 10: Complete control and evaluation record for the four-spin study. (a) The representative’s 96 physical control amplitudes. (b) GHZ fidelity and half-chain Rényi-2 entropy along its nominal evolution. (c) Fixed-pulse response to the specified coupling-error grid, with the median and observed range across all five initializations at each duration. All evaluated coupling offsets are shown.
Figure 11: Control-strength cost of reduced preparation error. Each line joins nominal and three-scenario designs from the same initialization. Across ten fixed confirmation points, the three-scenario designs have lower worst infidelity and higher control RMS in all three pairs. Lines identify paired initializations. All runs receive the same 660 single-scenario evaluation budget; Figure 4 b gives their individual error evaluations.
Grid
Design objective
J0
maxsJs
1922
Nominal
6.42×10−6
0.48828
1922
Five-scenario minimax
0.001209
0.46369
3842
Nominal
0.001578
0.51741
3842
Five-scenario minimax
0.004002
0.53757
Appendix
Table 9: Thermal-cloak objective trade-off. J0 is the nominal mean normalized squared error; maxsJs uses the five training material settings. Both designs remain frozen across grids. Lower is better.
Figure 12: Recovery of the 32-region FiPy conductivity field. True conductivity and the AD and FD recovered fields use the same dimensionless color scale and original 8×4 parameter zones. Relative L2 errors are 3.07×10−6 and 3.46×10−6 for AD and FD. These are the synthetic inverse-problem results, separate from thermal-cloak design.
Figure 13: Predictions at conditions excluded from fitting. (a) All 1,024 temperatures under FiPy’s withheld heater loading are compared with their synthetic targets. (b) Kwant LDOS profiles at energies −0.55 and 0.37 include all 36 sites, giving 72 withheld outputs. Both methods use the same model and discretization as the synthetic data. These checks test the fitted parameters at additional conditions, not robustness to model error.
Figure 14: Additional records for the two synthetic inverse problems. (a) FiPy objective relative to its initial value versus recorded optimizer time, using each method’s median-time process. The untimed initial value is omitted from the logarithmic time axis. (b) True and recovered onsite parameters at the 32 inferred sites of the Kwant LDOS problem. (c) All three Kwant optimizer times per method; vertical marks denote medians. Only total optimizer times, not iteration-level times, were recorded for Kwant. Held-out predictions appear in Figure 13 .
Differentiable programming makes programs and simulations amenable to gradient-based optimization by propagating derivatives through their computations.
Large language models (LLMs) have demonstrated abilities to perform complex tasks in multiple domains, including mathematical and scientific reasoning. In this work, we construct a novel agent framework for solving representative problems in scientific computing. The proposed agent incorporates a "rewriting-resolution-review-revision" logical chain via three reasoning LLMs (the Consultant, Reviewer, and Programmer) in a collaborative and interactive manner. The Consultant module enables knowledge transfer to link problems to professional domain insights, thereby rewriting problem descriptions through text augmentation. The Programmer module generates and executes well-structured code to solve problems. The Reviewer module enables self-debugging and self-refinement through interactive feedback with code runtime outputs. Through the end-to-end review mechanism, the Programmer's executable code undergoes iterative revision. We comprehensively evaluate the framework on typical partial differential equations (PDEs), engineering-oriented flow past a cylinder, ill-conditioned linear systems, and data-driven physical analysis problems. Compared with a single-model approach, this collaborative framework significantly improves the bug-free code generation rate and reduces non-physical solutions, establishing a reliable framework for autonomous code generation from natural language descriptions. The review mechanism improved the average execution success (bug-free code and non-NaN solutions) rate of DeepSeek R1 from 59% to 82%, while ChatGPT 4.1-mini improved from 66% to 87% and Gemini-2.5 improved from 60% to 84%. Automatic code generation and review provide a promising workflow for scientific computing; analysis of solver-development trajectories further reveals recurring implementation difficulties in boundary treatment, array indexing, and the preservation of numerical constraints.
Ao Cheng, Lei Zhang, Guowei He
The State Key Laboratory of Nonlinear Mechanics, Institute of Mechanics, Chinese Academy of Sciences, Beijing 100190, China · School of Engineering Science, University of Chinese Academy of Sciences, Beijing 100049, China
AI systems increasingly propose executable scientific models whose value depends on both their symbolic structure and their fitted continuous parameters. This makes parameter calibration the bottleneck of program-and-parameter co-search: an outer loop can generate thousands of candidate programs, but each needs an inner gradient-based optimization before it can be assessed. Staging each candidate into its own differentiable graph makes individual models fast but sacrifices the program-as-data property that keeps search fluid; interpreter-based approaches preserve programs as runtime data but pay interpreter overhead that dominates the numerical work. We present the Native Differentiable Virtual Machine (NDVM), a runtime representation that differentiates executable programs without compiling each candidate into a separate graph. NDVM separates symbolic structure from differentiable numeric state: tags, symbols, environments, and control remain native runtime data, while numeric payloads live in dense batched buffers with exact reverse-mode gradients recorded along the realized execution trace, so one evaluator walk is amortized across large populations of parameter vectors. A locked cost model of a real differentiable self-hosted Scheme interpreter motivates the design. We realize NDVM as a native runtime with forward and gradient equivalence to the reference backend, about 60x per-lane batch amortization, near-linear multicore scaling, and two independent front ends. In fixed-budget co-search over LLM-proposed programs, NDVM reaches high-quality solutions about 24x sooner in wall-clock time, suggesting runtime differentiation as a practical systems foundation for scientific discovery workflows.
Lucas Sheneman
Institute for Interdisciplinary Data Sciences University of Idaho
The boundary between program execution and gradient-based optimization has long limited the use of code itself as a learnable scientific model. We present a compiler that translates a self-hosting subset of Scheme into differentiable computation graphs for autograd backends. Because the subset can compile its own evaluator, this yields differentiable meta-circular interpretation (DMCI): a compiled Scheme interpreter executes programs supplied as data, while reverse-mode autodiff propagates gradients to continuous constants embedded in those programs. The interpreter is compiled once, so new programs inherit differentiability without recompilation or custom gradient machinery, while retaining closures, recursion, and data structures. We prove that gradients through the compiled interpreter are correct almost everywhere and show that they match direct compilation to numerical precision across 171 recursive and higher-order program-seed pairs. We then use DMCI for program-and-parameter co-search, where a large language model proposes Scheme programs and exact gradients calibrate their continuous parameters through a single frozen interpreter. This enables OpenEvolve-style program search in which an outer loop proposes discrete program structures and DMCI supplies exact gradient-based calibration of each candidate's continuous parameters. On battery capacity-fade data, the search recovers a knee-like degradation structure and improves held-out extrapolation over hand-crafted baselines on the harder early-extrapolation split, matching them on the later split. On a high-dimensional El Nino inverse problem, DMCI optimizes an interpreted Kalman-filter likelihood where gradient-free search fails. These results extend symbolic regression and neurosymbolic search from closed-form expressions to executable, stateful programs, making model-generated code directly optimizable against data.
Lucas Sheneman
Institute for Interdisciplinary Data Sciences University of Idaho