Performance: the step of the solvers' own difference quotients, and #1302's +4–36% wall on the fd path (Step 4 of #1002)

Step 4 of the plan on #1002, and the cost #1302 added, filed with #1002's close. Performance only: nothing here changes a result.

The step of the matrices that steer solvers (Step 4)

Where bngsim differences the state Jacobian to steer a solver (the steady-state march, the Newton polish, KINSOL), the sweep in include/bngsim/fd_jacobian.hpp steps each species by 1.5e-8 of the larger of itself and the state's largest concentration. Beside a species at 1e8, one near 0 is stepped by 1.5, and its column is a secant. In a matrix that only steers an iteration, that costs iterations, not correctness (the answers bngsim returns are confirmed or refused by #1100, #1294 and #1328). CVODE's own quotient floors each species' step by its error weight instead, which is the scale the solver cares about. Measuring whether that step saves iterations and wall time in the march and KINSOL on the corpus is the work.

What #1302 costs on the difference-quotient path

#1302 forms CVODE's dense quotient itself and tests each column for a jump: one more evaluation of f per tested column per Jacobian, three where a row may cross, one more where a column's mends are weighed. Steps and Jacobian counts are unchanged, and the closed-form default path is untouched. From #1302's measurement, jacobian="fd", codegen=False over the ODE perf suite, main → branch, interleaved best of 6 (evaluations / wall):

model	evaluations	wall
Repressilator	+21%	+4%
erk_activation	+32%	+10%
tcr_signaling	+36%	+13%
egfr_net_red	+39%	+18%
prion_aggregation (104 species)	+44%	+36%
SHP2_base_model (149, sparse)	+32%	+13%
fceri_ji (354, sparse)	+67%	+18%
egfr_net (356, sparse)	+78%	+27%

Round 1 of #1302's review suggested skipping the test for a model with no discontinuous operator, which would make the cost 0 for most models. It would not see a tanh steep enough to step across, so it was left. solver_stats["n_rhs_evals"] counts no Jacobian evaluation, so none of this shows there; counting them would make the cost visible to the nightly counters.

solver_stats: count a re-init's Jacobians once when no step follows it (

solver_stats: say what n_jac_formations counts, and hold the counts i…

solver_stats: count the RHS evaluations a difference-quotient Jacobia…

Merge pull request #1403 from Montana/perf/1333-fd-step-and-counts

Closing: neither part needs a further change.

Step 4 (the solvers' own step): measured on the corpus below. A different step saves no work in the default path, and in a forced Newton solve it loses more than it saves. The sweep's step stays.
Difference quotients that step across a jump in a rate law (#1002) #1302's cost: it is now visible. solver_stats: count the RHS evaluations a difference-quotient Jacobian costs (#1333) #1403 added solver_stats["n_jac_rhs_evals"] (the right-hand-side evaluations spent forming Jacobians) and n_jac_formations, and the nightly alerts on the first. The cost itself stays: Difference quotients that step across a jump in a rate law (#1002) #1302's review kept the jump test for the steep tanh it would otherwise miss.

Step 4: what was measured

An experimental build (not merged) replaced the matrices that steer the steady-state solvers where there is no closed form. On main those are KINSOL's own difference quotient, and the fd_jacobian.hpp sweep (1.5e-8·max(|y_j|, max_k |y_k|)) for the Newton distance that decides the polish. The build swapped in one of two alternatives:

ewt: a difference quotient stepped as CVODE steps its own, sqrt(eps)·max(|y_j|, atol_j);
dual: the exact Jacobian by dual numbers (A differenced Jacobian can be wrong at every step a check takes: take exact derivatives where there is no closed form (Step 3 of #1002) #1330's machinery), which has no step at all.

Corpus: the 585 ode_fullnet .net models and 1,323 rr_parity SBML models, jacobian="fd", models up to 400 species. All three variants finished on 1,432 models.

method="auto" (the default): every one of the 1,432 reaches its answer through the march. The march's Jacobian is CVODE's dense quotient, which already steps each species by its error weight. KINSOL steps by its own scale, not the largest concentration. The issue's premise, a step set by the largest species, holds only for the Newton distance. The alternatives change only what comes after it:

variant	steps	RHS evaluations	converged
ewt	+0.00%	+0.00%	+1 / −0
dual	+0.04%	+0.11%	+1 / −0

method="newton": 593 models settle by Newton with main's matrices.

variant	same method: steps	same method: RHS	lose Newton, fall back to the march	converged
ewt	−3.3% (577 models)	−3.0%	16 models (1–157 steps → 366–2,569)	+6 / −3
dual	−7.0% (581 models)	−6.4%	12 models (2–157 steps → 163–2,013)	+8 / −4

Where Newton still settles, the alternatives save a few percent of its iterations. They also lose Newton on 12–16 models where main's matrices found the root, and each of those costs hundreds of march steps. Summed over the corpus, both alternatives cost more steps than main's (dual +193%, ewt +300%, nearly all from those fallbacks). Wall time is a few milliseconds a model, and the variants ran in a fixed order, so it is not compared.

The sweep's step really is "a secant beside a species at 1e8". But a matrix that only steers an iteration finds the root as often with it as with an exact one, and on balance more often.

The experimental switch (BNGSIM_EXP_SOLVER_JAC=dual|ewt), on #1330's branch at c0fcb42; not merged

diff --git a/src/steady_state.cpp b/src/steady_state.cpp
index 99989b68..c6566cfc 100644
--- a/src/steady_state.cpp
+++ b/src/steady_state.cpp
@@ -105,6 +105,20 @@ static constexpr const char *kAutomaticDifferentiation = "automatic-differentiat
 class SteadyStateRhs {
   public:
     SteadyStateRhs(NetworkModel &model, const SteadyStateOptions &opts) : model_(model) {
+        // EXPERIMENT (issue #1333): the matrices that steer KINSOL and the
+        // Newton distance where there is no closed form.
+        if (const char *v = std::getenv("BNGSIM_EXP_SOLVER_JAC")) {
+            exp_mode_ = std::strcmp(v, "dual") == 0 ? 1 : std::strcmp(v, "ewt") == 0 ? 2 : 0;
+        }
+        {
+            const int ns0 = model.n_species();
+            exp_floor_.assign(static_cast<size_t>(ns0), opts.atol);
+            for (int i = 0; i < ns0; ++i) {
+                if (static_cast<size_t>(i) < opts.atol_vec.size()) {
+                    exp_floor_[static_cast<size_t>(i)] = opts.atol_vec[static_cast<size_t>(i)];
+                }
+            }
+        }
         const bool use_jit = !opts.codegen_c_source.empty();
         if (!use_jit && opts.codegen_so_path.empty()) {
             return; // interpreted
@@ -291,6 +305,38 @@ class SteadyStateRhs {
         return model_.fill_dual_dfdp(t, y, param_index, out, tie, tied);
     }
 
+    // EXPERIMENT (issue #1333).
+    int exp_mode() const { return exp_mode_; }
+    // The steering matrix with no closed form: dual numbers (mode 1), or a
+    // difference quotient stepped srur*max(|y_j|, atol_j) (mode 2). True where
+    // it filled `jac`.
+    bool fill_exp_solver_jacobian(double t, const double *y, double *jac) {
+        const int ns = model_.n_species();
+        if (exp_mode_ == 1 && fill_dual_jacobian(t, y, jac)) {
+            return true;
+        }
+        if (exp_mode_ == 0) {
+            return false;
+        }
+        std::vector<double> f0(static_cast<size_t>(ns)), f1(static_cast<size_t>(ns)),
+            yp(y, y + ns);
+        eval(t, y, f0.data());
+        for (int j = 0; j < ns; ++j) {
+            const double h_req =
+                1.4901161193847656e-8 * std::max(std::abs(y[j]), exp_floor_[static_cast<size_t>(j)]);
+            double xp = 0.0;
+            const double h = fd_probe(y[j], h_req, &xp);
+            yp[static_cast<size_t>(j)] = xp;
+            eval(t, yp.data(), f1.data());
+            yp[static_cast<size_t>(j)] = y[j];
+            for (int i = 0; i < ns; ++i) {
+                jac[static_cast<size_t>(j) * ns + i] = (f1[i] - f0[i]) / h;
+            }
+        }
+        return true;
+    }
+    const std::vector<double> &exp_floor() const { return exp_floor_; }
+
     // The analytical ∂f/∂p exists only in the compiled artifact — there is no
     // interpreted counterpart, so an absent symbol means the caller must
     // difference. Since GH #67 that is no longer "every Functional/MM model":
@@ -366,6 +412,8 @@ class SteadyStateRhs {
     CodegenUserDataForSO so_data_{};
     std::vector<double> param_buf_;
     std::string backend_ = "exprtk";
+    int exp_mode_ = 0;
+    std::vector<double> exp_floor_;
     // eval_dfdp scratch, kept on the object so a per-parameter loop does not
     // reallocate four n_species vectors per column.
     std::vector<double> zero_seed_, tmp1_, tmp2_, ydot_scratch_;
@@ -1798,7 +1846,12 @@ static int kinsol_reduced_jac_body(N_Vector y_ind, N_Vector /*fval*/, SUNMatrix
     // Every fill_dense_jacobian branch memsets its buffer, so resize (not assign)
     // is enough — the zeroing is not skipped, it is done once instead of twice.
     data->jac_full.resize(static_cast<size_t>(ns) * ns);
-    data->rhs->fill_dense_jacobian(0.0, data->jac_y_full.data(), data->jac_full.data());
+    if (data->rhs->has_analytical_jacobian()) {
+        data->rhs->fill_dense_jacobian(0.0, data->jac_y_full.data(), data->jac_full.data());
+    } else if (!data->rhs->fill_exp_solver_jacobian(0.0, data->jac_y_full.data(),
+                                                     data->jac_full.data())) {
+        return 1;
+    }
     ss_reduce_jacobian(data->jac_full.data(), ns, cl, cl.independent, data->jac_red);
     std::memcpy(SUNDenseMatrix_Data(J), data->jac_red.data(),
                 static_cast<size_t>(n_ind) * n_ind * sizeof(double));
@@ -1845,7 +1898,12 @@ static int kinsol_rhs(N_Vector y, N_Vector fval, void *ud) {
 static int kinsol_jac_body(N_Vector y, N_Vector /*fval*/, SUNMatrix J, void *ud, N_Vector /*tmp1*/,
                            N_Vector /*tmp2*/) {
     auto *data = static_cast<SteadyStateUserData *>(ud);
-    data->rhs->fill_dense_jacobian(0.0, N_VGetArrayPointer(y), SUNDenseMatrix_Data(J));
+    if (data->rhs->has_analytical_jacobian()) {
+        data->rhs->fill_dense_jacobian(0.0, N_VGetArrayPointer(y), SUNDenseMatrix_Data(J));
+    } else if (!data->rhs->fill_exp_solver_jacobian(0.0, N_VGetArrayPointer(y),
+                                                     SUNDenseMatrix_Data(J))) {
+        return 1;
+    }
     return 0;
 }
 
@@ -2370,7 +2428,9 @@ static SteadyStateResult solve_by_newton(NetworkModel &model, SteadyStateRhs &rh
     // single solve does not difference in one tier and use the closed form in
     // the other; the reduced callback carries the projection onto the unknown
     // set that the difference quotient it replaces got for free.
-    if (ss_install_solver_jacobian(rhs, opts)) {
+    if (ss_install_solver_jacobian(rhs, opts) ||
+        (rhs.exp_mode() != 0 && !(ss_want_analytical_jacobian(opts.jacobian) &&
+                                  rhs.has_analytical_jacobian()))) {
         flag = KINSetJacFn(kin_mem, use_reduced ? kinsol_reduced_jac : kinsol_jac);
         if (flag != KIN_SUCCESS) {
             throw std::runtime_error("KINSetJacFn failed");
@@ -3605,7 +3665,12 @@ static std::vector<double> ss_backward_columns(SteadyStateRhs &rhs, const double
     std::vector<double> y_pert(y, y + ns), f1(un), col(un);
     for (int j : cols) {
         const auto uj = static_cast<size_t>(j);
-        double step = fd_probe(y[j], -state_fd_step(y[j], y_scale), &y_pert[uj]);
+        double step = fd_probe(
+            y[j],
+            rhs.exp_mode() == 2
+                ? -1.4901161193847656e-8 * std::max(std::abs(y[j]), rhs.exp_floor()[uj])
+                : -state_fd_step(y[j], y_scale),
+            &y_pert[uj]);
         bool halved = false;
         if (y[j] >= 0.0 && y_pert[uj] < 0.0) {
             if (y[j] >= kSecantResolvable * std::abs(step)) {
@@ -3705,6 +3770,14 @@ static double ss_newton_distance(NetworkModel &model, SteadyStateRhs &rhs,
     const auto jacobian_at = [&](const double *at, std::vector<double> &Jm) {
         Jm.assign(un * un, 0.0);
         model.set_state_from(at);
+        if (!(want_analytical && rhs.has_analytical_jacobian()) && rhs.exp_mode() == 1 &&
+            rhs.fill_dual_jacobian(0.0, at, Jm.data())) {
+            return false; // exact
+        }
+        if (!(want_analytical && rhs.has_analytical_jacobian()) && rhs.exp_mode() == 2) {
+            rhs.fill_exp_solver_jacobian(0.0, at, Jm.data());
+            return true;
+        }
         return std::strcmp(ss_fill_state_jacobian(rhs, at, ns, sub, want_analytical, Jm.data(),
                                                   nullptr, /*exact=*/false),
                            "finite-difference") == 0;

The sweep ran steady_state() twice per variant on a fresh model each time, with jacobian="fd", and compared n_steps, n_rhs_evals, converged and method_used.

nightly: a flaky case's work is left out of the totals; list 385 and …

Navigation Menu

Repository navigation

Performance: the step of the solvers' own difference quotients, and #1302's +4–36% wall on the fd path (Step 4 of #1002) #1333

The step of the matrices that steer solvers (Step 4)

What #1302 costs on the difference-quotient path

Step 4: what was measured

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Performance: the step of the solvers' own difference quotients, and #1302's +4–36% wall on the fd path (Step 4 of #1002) #1333

Description

The step of the matrices that steer solvers (Step 4)

What #1302 costs on the difference-quotient path

Activity

wshlavacek commented on Nov 1, 2026

Step 4: what was measured

Metadata

Metadata

Assignees

Labels

Type

Projects

Milestone

Relationships

Development

Issue actions