Skip to content

fix: sampleCHRR returns NaN when the equality system is square - #696

Merged
edkerk merged 1 commit into
develop3from
fix/chrr-singular-equality-system
Aug 26, 2026
Merged

fix: sampleCHRR returns NaN when the equality system is square#696
edkerk merged 1 commit into
develop3from
fix/chrr-singular-equality-system

Conversation

@edkerk

@edkerk edkerk commented Aug 25, 2026

Copy link
Copy Markdown
Member

Problem

sampleCHRR fails on a model as small as the one shipped with the documentation
site (52 metabolites, 53 reactions):

sampleChebyshevCenter: LP infeasible - flux polytope has empty interior.

The polytope is not empty. The error comes from NaNs created much earlier.

sampleCHRR folds implicitly-fixed reactions into the equality system:

Aeq = [S; blk];          % blk: one row per reaction with ~zero FVA width
v0  = Aeq \ beq;         % particular solution

On this model exactly one reaction is fixed (ethIN, bounds [0 0]), so Aeq
is 53 x 53 — square — with rank 44. For a rank-deficient square system
MATLAB's backslash returns NaN, where for a rectangular one it would return a
basic solution. The NaNs flow into A_full/b_full, and the failure only
surfaces later, in sampleChebyshevCenter, as an infeasible LP.

The trigger is therefore nMets + nFixed == nRxns, which is not exotic — it is
one arithmetic coincidence away on any model.

Fix

v0 = lsqminnorm(Aeq, beq);

lsqminnorm returns the minimum-norm solution for a rank-deficient system of
any shape. As the surrounding comment already notes, v0 may be any particular
solution — the sampled distribution does not depend on the choice — so this
changes nothing except the degenerate case. It also matches the Python
implementation, which uses lstsq.

Verification

On the same model, before: the error above at every fixedWidthTol tried
(1e-7, 1e-6, 1e-4). After:

9 dimensions, 1 fixed, 53 x 200 samples, no NaNs
growth mean 0.0460, range 0.0097..0.0992

The Python implementation reports 9 dimensions and 1 fixed reaction for the
same model, so the two now agree on the geometry; the sampled growth
distribution matches too (0.0098..0.0911 there, a different RNG).

randomSampling's default achr path is untouched and still returns NaN-free
samples.

Found while writing the sampling page of the new user guide.

sampleCHRR builds Aeq = [S; rows fixing implicitly-determined reactions]
and takes a particular solution with `Aeq \ beq`. When the number of
metabolites plus the number of folded-in fixed reactions happens to equal
the number of reactions, Aeq is square; it is also rank-deficient, and
MATLAB's backslash returns NaN for that case rather than a least-squares
solution.

The NaNs propagate into the reduced polytope, and the failure surfaces
later and misleadingly as

    sampleChebyshevCenter: LP infeasible - flux polytope has empty interior.

Use lsqminnorm, which returns the minimum-norm solution for a
rank-deficient system of any shape -- and is what the Python reference
implementation already does, so the two now agree on the polytope.

Reproduces on the small yeast model shipped with the documentation site
(52 metabolites, 53 reactions, ethIN fixed at 0 -> Aeq is 53x53 with
rank 44). With the fix, sampling returns 9 dimensions and no NaNs,
matching the Python implementation's 9.
@github-actions

Copy link
Copy Markdown

Function test results

300 tests   274 ✅  1m 7s ⏱️
 25 suites   26 💤
  1 files      0 ❌

Results for commit ba36fc1.

@edkerk
edkerk merged commit 4ea9b09 into develop3 Aug 26, 2026
4 checks passed
@edkerk
edkerk deleted the fix/chrr-singular-equality-system branch August 26, 2026 07:19
edkerk added a commit to edkerk/raven-docs that referenced this pull request Aug 26, 2026
* Add user guide page 15: random sampling

Covers a first seeded sample, what a distribution says that a range does
not, the choice between interior sampling and random-objective vertices,
conditioning on a state before sampling, and the loopless screen behind
the objective list -- MATLAB and Python side by side.

Two things the examples measured that are worth calling out. The loop
pair from page 14 is sampled across nearly its whole 1000-unit range, so
loops are not a rare corner of the space but most of it; and sampled
growth averages under half the optimum, because near-uniform sampling
almost never lands on a vertex. Both make the same point: a sampled mean
describes the shape of the feasible space, not the organism.

The MATLAB CHRR block is skipped for now. randomSampling(...,'chrr')
fails on smallYeast with "flux polytope has empty interior" because
sampleCHRR's equality system is square and rank-deficient there, so its
particular solution is NaN -- SysBioChalmers/RAVEN#696. The documented
output is captured from a patched RAVEN; unskip once the fix ships.

* Print sampled statistics to a precision that survives the platform

CI put FRDS2's sampled minimum at 18.8 where this machine had 19.0, and
mean growth at 0.0497 against 0.0496. Same seed, same code: MATLAB's
samplers take their nullspace basis from null(), hence from LAPACK, so
the walk differs between Linux and Windows even though the distribution
does not. The Python chain reproduces exactly across both.

Print two decimals rather than four, and say all of this in 15.1 -- a
reader who seeds a chain and compares numbers with a colleague on
another machine should know which parts are the result.
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant