rafatahmed · Oct 20, 2026 · Oct 20, 2026 · Oct 20, 2026 · Oct 20, 2026
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -2,6 +2,8 @@
 
 ## Unreleased
 
+- Speed up the equilibrium solver and its searches (#37) with numerically equivalent results (largest change across all PHREEQC fixtures: 6e-12 relative in species, 4e-13 in the nitric dose). The Jacobian and component sums now use each species' non-zero stoichiometry instead of every component (94% of runtime). The nitric-target, acid-option and alkalinity-to-carbon searches share one bracketed Illinois root finder with the same stopping rules as before, in place of fixed bisections. The calculated-pH search keeps its bisection, so the shipped phosphate conductivity fit (#37) reproduces exactly. The mixed-product nitric target fell from 44 solves and about 70 s to 12 solves and under 2 s (median 1.6 s) on the recording laptop. Benchmark v1 adds `track_search_solves`, which reports the solves per call.
+- Fix the equilibrium solver raising 'line search stalled' at some valid states: after a cold start a near-singular Jacobian (seen with Fe-EDDHA in the mixed-product solution at pH 7.396) gave an unusable Newton step. A damped least-squares (Levenberg-Marquardt) step now recovers. It runs only where the solver previously failed, so converging solves are unchanged.
 - Add Potassium Chloride (KCl) and Calcium Chloride (Dihydrate) (CaCl2.2H2O) to the catalogue (30 products), with identity-theoretical assays (K 52.4469 / Cl 47.5531 %; Ca 27.2626 / Cl 48.2289 %) sourced to Sigma-Aldrich P3911 and 223506. Calcium chloride is the first calcium source without nitrate (#39). Chloride needs consent: neither product is used unless the caller targets Cl, allows it (`allow_ions`) or sets a Cl maximum, so no default recipe changes.
 - Record the ChemicalBook technical data sheets (revision 2026-04-16, SHA-256 recorded) as sources for Ammonium Sulfate and Urea; they confirm identity and purity (laboratory grade). Assays are unchanged.
 - `recommend(..., maximum_concentrations={...})` and `plan_fertilizer_workflow` now use incidental-ion ceilings in the fit: an ion with a positive ceiling (for example `{"Cl": 70}`) is admitted, is not penalised up to the ceiling, and is pulled back when the fit exceeds it. Any remaining excess is reported in `limitViolations` with `feasible` false. Without ceilings the result is unchanged.

diff --git a/benchmarks/performance.py b/benchmarks/performance.py
@@ -121,6 +121,23 @@ def peakmem_solve(self, fixture):
     def track_final_state_iterations(self, fixture):
         return self.result.iterations
 
+    def track_search_solves(self, fixture):
+        # Total equilibrium solves per call; an acid target runs a search, so its
+        # final-state iterations alone understate the work.
+        from flahax import aqueous_model
+        original, count = aqueous_model.solve, [0]
+
+        def counted(*args, **kwargs):
+            count[0] += 1
+            return original(*args, **kwargs)
+
+        aqueous_model.solve = counted
+        try:
+            self.solve(self.data)
+        finally:
+            aqueous_model.solve = original
+        return count[0]
+
     def teardown(self, fixture):
         require_chemistry(self.data, self.result, self.acid)
 

diff --git a/docs/site/evidence/benchmark-v1.md b/docs/site/evidence/benchmark-v1.md
@@ -67,7 +67,9 @@ absolute/relative molality and acid-dose tolerances, absolute log-activity and
 saturation-index tolerances, and phase-amount tolerances. Conservation uses
 `1e-15 + 1e-9 * component total`; phase complementarity uses `1e-7` SI.
 Report every comparison and available final-state iteration count, including
-failed comparisons. Do not silently widen tolerances to improve results.
+failed comparisons. `track_search_solves` counts every equilibrium solve per call,
+so an acid target's search work is reported separately from the iterations of
+its final state. Do not silently widen tolerances to improve results.
 
 These captures prescribe pH. The nitric case solves reagent demand to a
 prescribed target; it does not independently validate prediction of an unknown
@@ -151,6 +153,9 @@ timed samples. The equilibrium benchmark has a 900-second process timeout;
 this budget includes setup. A timeout is a benchmark-budget failure, not proof
 that one solve exceeded 900 seconds or failed to converge. Preserve it as a
 missing measurement and use a separately identified run if changing the budget.
+The 2026-10-20 snapshots predate the solver speed-up for #37 (sparse Jacobian and
+bracketed root finding), which cut the acid-target call from about 70 s to under
+2 s on the recording laptop; compare only runs from the same revision.
 
 The manual `Benchmark v1` CI workflow retains JSON, ASV records, and its
 measurement log as artifacts. A correctness failure makes the workflow fail.

diff --git a/src/flahax/acid_selection.py b/src/flahax/acid_selection.py
@@ -12,6 +12,7 @@
 from .delivery_contracts import DeliveryError
 from .delivery_quantities import rescore_with_reagent, MilligramsPerLitre
 from .nutrient_acceptance import assess_incidental
+from .root_finding import shrink_bracket
 
 # Pure, anhydrous acid identities; commercial aqueous products supply their own
 # mass fraction/density. Hydrated citric solids require a separately defined stock.
@@ -119,18 +120,18 @@ def at(dose):
                 if target.charge_balance < initial.charge_balance - 1e-12:
                     raise DeliveryError('wrong_reagent_direction', 'target requires base under this phase policy')
                 low, high = 0., 1e-6
+                g_low = target.charge_balance - initial.charge_balance
                 for _ in range(32):
-                    if at(high).charge_balance <= initial.charge_balance:
+                    g_high = at(high).charge_balance - initial.charge_balance
+                    if g_high <= 0:
                         break
+                    low, g_low = high, g_high
                     high *= 2
                 else:
                     raise DeliveryError('nonconvergent', 'acid demand could not be bracketed')
-                for _ in range(40):
-                    middle = (low + high) / 2
-                    if at(middle).charge_balance > initial.charge_balance:
-                        low = middle
-                    else:
-                        high = middle
+                # Same final width as 40 bisection steps from [0, high].
+                low, high = shrink_bracket(lambda d: at(d).charge_balance - initial.charge_balance,
+                                           low, g_low, high, g_high, width=high * 2. ** -40)
                 dose = (low + high) / 2
                 target = at(dose)
             grams = dose * mass * molar_mass / fraction

diff --git a/src/flahax/aqueous_model.py b/src/flahax/aqueous_model.py
@@ -11,6 +11,7 @@
 import json
 import math
 from .delivery_contracts import DeliveryError
+from .root_finding import shrink_bracket
 
 LN10 = math.log(10.)
 
@@ -50,6 +51,28 @@ def linear(a, b):
         x[i]=(a[i][n]-sum(a[i][j]*x[j] for j in range(i+1,n)))/a[i][i]
     return x
 
+def _damped_step(evaluate, x, residual, matrix):
+    """Levenberg-Marquardt step that lowers the squared residual, or None."""
+    m=len(residual)
+    normal=[[sum(matrix[r][i]*matrix[r][j] for r in range(m)) for j in range(m)] for i in range(m)]
+    gradient=[sum(matrix[r][i]*residual[r] for r in range(m)) for i in range(m)]
+    scale=[max(normal[i][i],1e-300) for i in range(m)]
+    merit=sum(v*v for v in residual)
+    damping=1e-3
+    for _ in range(40):
+        system=[[normal[i][j]+(damping*scale[i] if i==j else 0.) for j in range(m)] for i in range(m)]
+        try:
+            step=linear(system,[-g for g in gradient])
+        except DeliveryError:
+            step=None
+        if step is not None:
+            trial=[v+d for v,d in zip(x,step)]
+            got=evaluate(trial)
+            if got and sum(v*v for v in got[0])<merit:
+                return trial
+        damping*=10
+    return None
+
 def solve(totals, ph, *, phases=(), ionic_strength=None, temperature_c=25.):
     if temperature_c != 25:
         raise DeliveryError('temperature_out_of_range','reaction dataset is restricted to 25 C')
@@ -71,6 +94,8 @@ def solve(totals, ph, *, phases=(), ionic_strength=None, temperature_c=25.):
         raise DeliveryError('unknown_phase','phase not in the pinned reaction set')
     candidates={p:all_phases[p] for p in phases if set(all_phases[p]['powers'])<=permitted}
     stoich=[[r['powers'].get(b,0.) for b in bases] for r in records]
+    # Each species involves few components; iterate only its non-zero stoichiometry.
+    sparse=[[(j,v) for j,v in enumerate(s) if v] for s in stoich]
     strength=ionic_strength if ionic_strength is not None else min(.05,max(1e-8,sum(totals.values())))
     logs=[math.log10(totals[b])-3 for b in bases]
     active=[]
@@ -84,8 +109,14 @@ def lgamma(z):
     for outer in range(100):
         offsets=[r['beta']-ph*r['powers'].get('H+',0)+math.log10(water)*r['powers'].get('H2O',0)-lgamma(r['charge']) for r in records]
         def evaluate(x, jacobian=False):
-            amounts=[10.**max(-300,min(100,o+sum(a*v for a,v in zip(s,x[:n])))) for o,s in zip(offsets,stoich)]
-            sums=[sum(a*s[j] for a,s in zip(amounts,stoich))+sum(x[n+k]*all_phases[p]['powers'].get(b,0.) for k,p in enumerate(active)) for j,b in enumerate(bases)]
+            amounts=[10.**max(-300,min(100,o+sum(v*x[j] for j,v in nz))) for o,nz in zip(offsets,sparse)]
+                sums=[0.]*n
+            for a,nz in zip(amounts,sparse):
+                for j,v in nz:
+                    sums[j]+=a*v
+            for k,p in enumerate(active):
+                for j,b in enumerate(bases):
+                    sums[j]+=x[n+k]*all_phases[p]['powers'].get(b,0.)
             if any(v<=0 for v in sums):
                 return None
             residual=[math.log10(v/totals[b]) for b,v in zip(bases,sums)]
@@ -94,9 +125,16 @@ def evaluate(x, jacobian=False):
                 residual.append(sum(r['powers'].get(b,0)*x[j] for j,b in enumerate(bases))-ph*r['powers'].get('H+',0)+math.log10(water)*r['powers'].get('H2O',0)-r['beta'])
             if not jacobian:
                 return residual,amounts
+            moments=[[0.]*n for _ in range(n)]
+            for a,nz in zip(amounts,sparse):
+                for j,vj in nz:
+                    weighted=a*vj
+                    row=moments[j]
+                    for k,vk in nz:
+                        row[k]+=weighted*vk
             matrix=[]
             for j,b in enumerate(bases):
-                row=[sum(a*s[j]*s[k] for a,s in zip(amounts,stoich))/sums[j] for k in range(n)]
+                row=[value/sums[j] for value in moments[j]]
                 row += [all_phases[p]['powers'].get(b,0.)/(sums[j]*LN10) for p in active]
                 matrix.append(row)
             for p in active:
@@ -120,7 +158,12 @@ def evaluate(x, jacobian=False):
                         break
                     alpha*=.5
                 else:
-                    raise DeliveryError('nonconvergent',f'equilibrium line search stalled (residual {norm:g})')
+                    # A near-singular Jacobian can give a useless Newton direction
+                    # (seen after a cold start with a strong chelate); take a damped
+                    # least-squares step instead. Paths that converge today are unchanged.
+                    x=_damped_step(evaluate,x,residual,matrix)
+                    if x is None:
+                        raise DeliveryError('nonconvergent',f'equilibrium line search stalled (residual {norm:g})')
             else:
                 raise DeliveryError('nonconvergent','aqueous Newton iteration limit')
             logs=x[:n]
@@ -179,14 +222,12 @@ def at(dose):
     target=at(0.)
     if target.charge_balance < initial.charge_balance-1e-12:
         raise DeliveryError('wrong_reagent_direction','target requires base')
-    while at(hi).charge_balance>initial.charge_balance:
+    f_lo=target.charge_balance-initial.charge_balance
+    while (f_hi:=at(hi).charge_balance-initial.charge_balance)>0:
+        lo,f_lo=hi,f_hi
         hi*=2
-    for _ in range(40):
-        mid=(lo+hi)/2
-        target=at(mid)
-        if target.charge_balance>initial.charge_balance:
-            lo=mid
-        else:
-            hi=mid
+    # Same final width as 40 bisection steps from [0, hi], in far fewer solves.
+    lo,hi=shrink_bracket(lambda dose:at(dose).charge_balance-initial.charge_balance,
+                         lo,f_lo,hi,f_hi,width=hi*2.**-40)
     dose=(lo+hi)/2
     return AcidResult(initial,at(dose),dose)
diff --git a/src/flahax/measured_water.py b/src/flahax/measured_water.py
@@ -14,6 +14,7 @@
 
 from .aqueous_model import AqueousResult, chemistry, solve
 from .delivery_contracts import DeliveryError
+from .root_finding import shrink_bracket
 
 BALANCE_IONS = {'Cl-': -1, 'Na+': 1}
 DEFAULT_MAXIMUM_IMBALANCE_PERCENT = 5.
@@ -75,19 +76,14 @@ def alk(carbon):
     if alkalinity <= base:
         raise DeliveryError('alkalinity_out_of_range',
                             'measured alkalinity is not above the carbon-free alkalinity at the measured pH')
-    low, high = 0., 1e-6
-    while alk(high) < alkalinity:
+    low, g_low, high = 0., base - alkalinity, 1e-6
+    while (g_high := alk(high) - alkalinity) < 0:
+        low, g_low = high, g_high
         high *= 4
         if high > .1:
             raise DeliveryError('alkalinity_out_of_range', 'alkalinity requires more carbon than the Davies domain allows')
-    for _ in range(200):
-        middle = (low + high) / 2
-        if alk(middle) < alkalinity:
-            low = middle
-        else:
-            high = middle
-        if high - low <= 1e-12 * high:
-            break
+    low, high = shrink_bracket(lambda carbon: alk(carbon) - alkalinity, low, g_low, high, g_high,
+                               relative=1e-12)
     return (low + high) / 2
 
 

diff --git a/src/flahax/root_finding.py b/src/flahax/root_finding.py
@@ -0,0 +1,47 @@
+"""Bracketed root finding shared by the acid-dose and alkalinity searches.
+
+Each search previously bisected a fixed number of times, solving the full
+equilibrium model at every step. Illinois false position keeps the same sign
+bracket and the same stopping rule but converges superlinearly on these smooth
+monotonic functions. A step outside the bracket, or two steps that fail to halve
+it, fall back to bisection, so the worst case stays within twice bisection's.
+"""
+
+
+def shrink_bracket(g, lo, g_lo, hi, g_hi, *, width=0., relative=0., accept=None, max_steps=200):
+    """Narrow [lo, hi], where g changes sign, until hi - lo <= width + relative * hi.
+
+    ``g_lo`` and ``g_hi`` are g(lo) and g(hi) with opposite signs (either may be
+    zero). ``accept(x, g_x, lo, hi)`` may end the search early at an evaluated
+    point. Returns the final (lo, hi); both equal x when g(x) is exactly zero.
+    """
+    if g_lo == 0:
+        return lo, lo
+    if g_hi == 0:
+        return hi, hi
+    side = 0
+    widths = [hi - lo]
+    for _ in range(max_steps):
+        if hi - lo <= width + relative * abs(hi):
+            break
+        x = hi - g_hi * (hi - lo) / (g_hi - g_lo) if g_hi != g_lo else (lo + hi) / 2
+        # Guarantee progress: bisect when two steps have not halved the bracket.
+        if not lo < x < hi or (len(widths) >= 3 and widths[-1] > widths[-3] / 2):
+            x = (lo + hi) / 2
+        g_x = g(x)
+        if g_x == 0:
+            return x, x
+        if (g_x > 0) == (g_lo > 0):
+            lo, g_lo = x, g_x
+            if side == -1:
+                g_hi /= 2
+            side = -1
+        else:
+            hi, g_hi = x, g_x
+            if side == 1:
+                g_lo /= 2
+            side = 1
+        widths.append(hi - lo)
+        if accept is not None and accept(x, g_x, lo, hi):
+            return x, x
+    return lo, hi
diff --git a/tests/test_chemistry_acceptance.py b/tests/test_chemistry_acceptance.py
@@ -29,6 +29,17 @@ def test_urea_and_potassium_carbonate_stoichiometry(self):
         carbonate=catalogue_dose_totals({'Potassium Carbonate':1.})
         self.assertEqual(carbonate['K+'],2*carbonate['CO3-2'])
 
+    def test_near_singular_newton_step_does_not_stall(self):
+        # Cold start at this pH met a near-singular Jacobian (Fe-EDDHA) and raised
+        # 'line search stalled'; a damped least-squares step now recovers.
+        from benchmarks.cases import read_fixture
+        totals=read_fixture('products/mixed_target_ph_nitric')['runtime_totals']
+        for ph in (7.3962594096730125,7.406259409673012):
+            with self.subTest(ph=ph):
+                result=solve(totals,ph)
+                for component,residual in result.residuals.items():
+                    self.assertLess(abs(residual),1e-9*result.totals[component]+1e-15)
+
     def test_chloride_salts_keep_chloride_and_add_no_nitrate(self):
         calcium=catalogue_dose_totals({'Calcium Chloride (Dihydrate)':1.})
         self.assertNotIn('NO3-',calcium)

diff --git a/tests/test_root_finding.py b/tests/test_root_finding.py
@@ -0,0 +1,48 @@
+"""Shared bracketed root finding used by the dose, pH and alkalinity searches."""
+import unittest
+
+from flahax.root_finding import shrink_bracket
+
+
+class ShrinkBracketTests(unittest.TestCase):
+    def test_reaches_the_requested_width_in_few_evaluations(self):
+        calls = []
+
+        def g(x):
+            calls.append(x)
+            return 2. - x * x
+
+        lo, hi = shrink_bracket(g, 0., 2., 2., -2., width=2. * 2. ** -40)
+        self.assertLessEqual(hi - lo, 2. * 2. ** -40)
+        self.assertAlmostEqual((lo + hi) / 2, 2 ** .5, places=12)
+        self.assertLess(len(calls), 20)  # bisection needs 40
+
+    def test_relative_width_and_decreasing_functions(self):
+        lo, hi = shrink_bracket(lambda x: 1e-3 - x, 0., 1e-3, 1., 1e-3 - 1., relative=1e-12)
+        self.assertLessEqual(hi - lo, 1e-12 * hi + 1e-18)
+        self.assertAlmostEqual(lo, 1e-3, delta=1e-14)
+
+    def test_exact_zero_and_zero_endpoint(self):
+        self.assertEqual(shrink_bracket(lambda x: x - .5, 0., -.5, 1., .5, width=1e-9), (.5, .5))
+        self.assertEqual(shrink_bracket(lambda x: x, 0., 0., 1., 1.), (0., 0.))
+
+    def test_accept_ends_at_an_evaluated_point(self):
+        seen = []
+
+        def g(x):
+            seen.append(x)
+            return x - .3
+
+        lo, hi = shrink_bracket(g, 0., -.3, 1., .7, accept=lambda x, gx, lo, hi: abs(gx) < 1e-6)
+        self.assertEqual(lo, hi)
+        self.assertIn(lo, seen)
+
+    def test_stray_secant_steps_fall_back_to_bisection(self):
+        # A steep step function defeats false position; the bracket must still shrink.
+        lo, hi = shrink_bracket(lambda x: -1. if x < .123456 else 1e9, 0., -1., 1., 1e9, width=1e-9)
+        self.assertLessEqual(hi - lo, 1e-9)
+        self.assertTrue(lo <= .123456 <= hi)
+
+
+if __name__ == "__main__":
+    unittest.main()