cs.AIOct 3, 2026
SaveRecursive Improvement of a Differentiable Scientific Software Ecosystem
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 ·
Abstract
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. | 14 composed cases pass reference-value and finite-difference checks (Appendix C.4 ). |
| Transport and flux sensitivities OpenPNM | Pressure-driven flow, diffusion, heat, reaction–diffusion and advection–diffusion sensitivities. | 66/80 pass value/derivative checks; 14 hard-threshold invasion cases skipped (Appendix C.4 ). |
| Simulator-connected sensitivities ParFlow | Pressure and saturation sensitivities; groundwater-head and water-storage objective gradients. | 80 sensitivity cases pass original-simulator finite-difference checks (Appendix B.1 ). |
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 : FD evaluates derivatives, while AD returns a gradient over all independent parameters. Solid bars have ; hatched bars have . Bars show ratios of median runtimes, , 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 and encoded phase . (b) Error for three paired initializations, where 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 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 grid and evaluated at 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 |
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 |
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 |
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 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 | with and a controlled 12-dimensional Hamiltonian; 144 real Hermitian-basis controls. | All 144 active controls; base and active set vary across the smaller- series. |
| QuSpin | Real target-state overlap after symmetric-Hamiltonian evolution; 528 independent real matrix coefficients. | First 256 symmetric directions. |
| BioDynaMo | Weighted output of one periodic diffusion step; 4,096 input concentrations in an extension kernel. | Seeded 256 voxels. |
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 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 | 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 |
Table 6: Definitions and reported timing medians for all workloads in Figure 2 . gives derivatives requested by FD / the complete task-level gradient returned by AD (after contraction for MCell); when , FD requests a selected subset contained in the full AD result, not the same complete gradient. is the source-reported repeat count, not independent processes. Times are milliseconds and follow the Figure 2 ratio order. MCell ( ) uses five requests for timing medians; two 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 mesh; Kwant grows to 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, , 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 |
| 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: 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 | Package rule differentiating Hamiltonian-active evolution | Eight programs migrated; two more pass ( /47); later control reuses the rule |
| QuantumOptics: compose further inference and control tasks | Nine programs using existing evolution, steady and spectral rules | Workflow retests: of 49 pass; no new kernel implied |
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 |
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 | ||
|---|---|---|---|
| Nominal | 0.48828 | ||
| Five-scenario minimax | 0.001209 | 0.46369 | |
| Nominal | 0.001578 | 0.51741 | |
| Five-scenario minimax | 0.004002 | 0.53757 |
Table 9: Thermal-cloak objective trade-off. is the nominal mean normalized squared error; 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 parameter zones. Relative errors are and 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 and 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 .
More in this topic
Differentiable Programming
Differentiable programming makes programs and simulations amenable to gradient-based optimization by propagating derivatives through their computations.
See all 45 papers