Solve binding ceilings as constraints; test NNLS optimality per column (#68) by rafatahmed · Pull Request #69 · rafatahmed/Flahax

rafatahmed · 2026-10-20T23:25:20Z

Closes #68.

Problem

With the default ammonium limit (#54), recommend() was 2–5× slower than without it, even where the limit did not bind. On the Benchmark v1 formulation workload it was 31.8 ms against 14.5 ms before #54 (same machine). Investigating it turned up a correctness issue as well: recipes under a binding ceiling depended on the solver's starting point.

Root causes (each measured against SciPy NNLS on every fit)

NNLS stopped early. The optimality test (from Investigate NNLS objective discrepancy for near-dependent columns #36) used one tolerance for every column, 10·max(m,n)·ε·max|A|·‖b‖. A heavily weighted ceiling row inflated it for every column, including columns that row does not touch. On the 585 ceiling fits of the benchmark, FlahaX left products at zero whose dual weight was still positive (up to 1.9e-6). Doses were up to 4.2e-3 g/L away from the optimum.
A binding ceiling was a 1024-weight penalty. That problem is so flat in double precision that its doses are fixed only to about 2e-3 g/L, for SciPy too. This is also why its fits fell back to Householder: 15 fallbacks per call on the benchmark target.

Fix (pure Python)

NNLS: each column's dual weight cross[j] − Σ G[j][k]·x[k] is compared with the rounding error of its own terms, 10·max(m,n)·ε·(Σ|a_rj·b_r| + Σ|G[j][k]|·x[k]). This keeps Investigate NNLS objective discrepancy for near-dependent columns #36's scale invariance (tested at scales 1e-8, 1, 1e8). The bound comes out of the same loop as the weights.
Ceilings as constraints:
- A binding ceiling gets a weight-8 row aimed at aim − shift. The shift is the scaled Lagrange multiplier, and the aim is the limit less 1 ppm.
- The ion's excess is non-increasing in the shift, and affine while NNLS keeps its active set. So the shift is found by a bracketed search (secant extrapolation, then Illinois, with a bisection safeguard), to 1e-9 of the limit.
- A ceiling whose multiplier falls to zero is released. A ceiling that cannot be met, such as water alone above the limit, is still reported in limitViolations.
Gram reuse: the matrix is the same throughout a ceiling search, so nnls(..., gram=) shares one Gram cache.

Evidence

Check	`main` (`4dffdf8`)	This PR
Ceiling fits vs SciPy, worst dose difference	4.2e-3 g/L	1.4e-10 g/L
Warm vs cold start, worst dose difference	4.2e-3 g/L	6.6e-9 g/L
Householder fallbacks (104 runs)	2,004	255
52 recipes, default ammonium limit	4.20 s	2.78 s (−34%)
52 recipes, Cl ceiling 3 mmol/L	3.77 s	3.87 s (+3%)
ASV formulation workload, Ca 0 / 20 (same laptop)	30 / 17 ms	16 / 13 ms
Complete recipes, default	34/52	34/52
Complete recipes, Cl ceilin 8000 g	41/52	42/52
Targets passing, default / Cl ceiling	571 / 598	571 / 599
Limit violations	0	0

Changed recipes: 10 of the 104 runs change (5 with the default, 5 with the Cl ceiling). Each has an equal or better target fit; the worst error goes from 0.004 % to 0.002 %, or from 0.002 % to 0.000 %. The new recipes are the unique constrained optimum. The old ones were arbitrary points inside the flat band.
Benchmark v1 on this branch: run 39514737671 succeeded. Formulation measured 11.6 / 9.5 ms; that's a different runner from the published snapshot, so compare the same-laptop row above. Locally, NNLS is 7/7 and chemistry 21/21.

Tests

371 tests pass (1 skipped).
tests/test_ceiling_constraint.py:
- a small NNLS problem where a large row hid a small column (on main it returns x₂ = 0 instead of 1);
- scale invariance;
- Generic Grow's binding ammonium limit lands on its aim (fails on main);
- warm and cold starts agree.
The usage guide and the recommend() docstring now describe ceilings as constraints. The API reference was regenerated.

…lumn (#68) recommend() with the default ammonium limit was 2-5x slower than without it, even where the limit did not bind, and its recipes depended on the solver's starting point. Measured against SciPy on every fit: - NNLS stopped early. Its optimality test used one tolerance scaled by the whole matrix (max|A| |b|); a heavily weighted ceiling row inflated it for every column, hiding real descent directions (a zero product with a positive dual weight of 1.9e-6; doses off by up to 4.2e-3 g/L). Each column's dual weight is now compared with the rounding error of its own terms. This keeps the scale invariance #36 needed. - A binding ceiling was a 1024-weight penalty. That problem is so flat in double precision that doses are fixed only to ~2e-3 g/L, for SciPy too. The ceiling is now a constraint: a weight-8 row aimed at aim - shift, where the shift (the scaled Lagrange multiplier) is found by a bracketed secant/Illinois search; the excess is non-increasing in the shift and affine while NNLS keeps its active set. The matrix is the same across the search, so one Gram cache is shared. Every ceiling fit now matches SciPy to 1.4e-10 g/L and warm and cold starts agree to 7e-9 g/L. Over the 52 benchmark recipes recommend() is 34% faster with the default ammonium limit (4.20 -> 2.78 s) and unchanged with a 3 mmol/L Cl ceiling (3.77 -> 3.87 s); Benchmark v1's formulation workload drops from about 30 to 16 ms. Default pass counts are unchanged (34/52, 571/636); with the Cl ceiling 41 -> 42 recipes and 598 -> 599 targets; no limit violations. Benchmark v1 scorecards: NNLS 7/7, chemistry 21/21.

read-the-docs-community · 2026-10-20T23:26:28Z

Documentation build overview

📚 Flahax | 🛠️ Build #35283688 | 📁 Comparing 16cbee0 against latest (4dffdf8)

🔍 Preview build

3 files changed

± changelog.html
± api/reference.html
± guide/usage.html

Hosted run 39515183621 on main at df0d30b (Benchmark v1, success; normal CI run 39515081885, success). Every correctness gate passes: 7/7 NNLS, 21/21 stored chemistry, 17/17 timing workloads; complete recipes 34/52. The unmerged 4dffdf8 snapshot predates #69 and is removed; its capture notes record why.

Solve binding ceilings as constraints; test NNLS optimality per column (#68)

Hosted run 39515183621 on main at df0d30b (Benchmark v1, success; normal CI run 39515081885, success). Every correctness gate passes: 7/7 NNLS, 21/21 stored chemistry, 17/17 timing workloads; complete recipes 34/52. The unmerged 4dffdf8 snapshot predates #69 and is removed; its capture notes record why.

rafatahmed merged commit df0d30b into main Oct 20, 2026
7 checks passed

rafatahmed deleted the fix/ammonium-ceiling-speed branch October 20, 2026 23:26

This was referenced Oct 20, 2026

Publish Benchmark v1 snapshot df0d30b: every correctness gate passing (#66) #67

Merged

test: benchmark current solver convergence and numerical performance #42

Closed

rafatahmed self-assigned this Oct 23, 2026

rafatahmed added a commit that referenced this pull request Nov 7, 2026

Merge pull request #69 from rafatahmed/fix/ammonium-ceiling-speed

c9b2c53

Solve binding ceilings as constraints; test NNLS optimality per column (#68)

Navigation Menu

Repository navigation

Solve binding ceilings as constraints; test NNLS optimality per column (#68) - #69

Solve binding ceilings as constraints; test NNLS optimality per column (#68)#69
rafatahmed merged 1 commit into
mainfrom
fix/ammonium-ceiling-speed

Uh oh!

Uh oh!

Uh oh!

Reviewers

Uh oh!

Assignees

Labels

Projects

Milestone

Development

Uh oh!

1 participant

Conversation

Problem

Root causes (each measured against SciPy NNLS on every fit)

Fix (pure Python)

Evidence

Tests

Uh oh!

Documentation build overview

Uh oh!

Uh oh!

Reviewers

Uh oh!

Assignees

Labels

Projects

Milestone

Development

Uh oh!

1 participant