`SteadyStateProblem(FSPSystem, ...)` is underdetermined and `SSRootfind` can return the zero vector

Summary

The finite-state steady-state formulation currently solves the homogeneous system

Q u = 0.

This system admits both the stationary distribution and the zero vector. SSRootfind() can therefore return a near-zero vector with ReturnCode.Success; the test then divides that vector by its near-zero sum, amplifying roundoff into a non-stationary signed vector.

I reproduced this on an exact, unmodified checkout of current main (5c031a11ba0243d992fb6210745e598c0be37f2a) with a clean, dedicated depot. This is not caused by the generated RHS or a permutation/matrix construction error: the generated steady RHS agrees exactly with the sparse matrix, and the permuted matrix is exactly similar to the original matrix in the deterministic probe.

Clean-main reproduction

Environment:

• Julia 1.12.6
• FiniteStateProjection 0.4.0 at 5c031a11ba0243d992fb6210745e598c0be37f2a
• Catalyst 16.2.2
• ModelingToolkitBase 1.55.0
• SciMLBase 3.36.1
• OrdinaryDiffEq 7.1.3
• SteadyStateDiffEq 2.11.3

Command:

GROUP=Core JULIA_PKG_PRECOMPILE_AUTO=0 julia +1.12 --startup-file=no --history-file=no --color=no --project=. -e 'using Pkg; Pkg.test()'

Result:

Test Summary:        | Pass  Fail  Total      Time
Core/birthdeath2D.jl |   15     3     18  18m57.7s
RNG of the outermost testset: Random.Xoshiro(0xdeee60cefff28b4e, 0x5854311b29c3d09f, 0xe2f99810f89bbe6c, 0xec3548543c05b3ba, 0xcc3fe6c2fc709f5d)

The failures are the stationary marginal assertions at test/birthdeath2D.jl:74-75 and the permuted stationary solution comparison at line 116.

Deterministic exact replay

I also replayed a second printed failing RNG state exactly:

Random.Xoshiro(
    0xda3d877b33ec4332,
    0x27a3fa717f3046c6,
    0x79f6bd931f3114ea,
    0x14f421f56ab82713,
    0x45a81865387b919b,
)

Its first four draws produce this parameter map:

pmap = [
    :r1 => 3.142200691084267,
    :r2 => 0.494653851350257,
    :s1 => 3.4491822160853136,
    :s2 => 1.5870402869698628,
]

With Nmax = 45, the exact diagnostics are:

transient_rhs_matrix_max_abs_diff=0.0
steady_rhs_matrix_max_abs_diff=0.0
steady_matrix_size=(2116, 2116)
steady_matrix_nnz=10396
steady_matrix_column_sum_inf=1.1546319456101628e-14

fsp_run_1.retcode=Success
fsp_run_1.raw_sum=-2.3337091563614036e-16
fsp_run_1.raw_norm_inf=1.7505893086657684e-17
fsp_run_1.raw_residual_inf=1.1596592831945752e-16
fsp_run_1.normalized_min=-0.057677175833832924
fsp_run_1.normalized_max=0.07501317393788672
fsp_run_1.normalized_residual_inf=0.49691679875081557
fsp_run_1.marginal_a_max_abs_error=0.12911011058932598
fsp_run_1.marginal_b_max_abs_error=0.12767355428800362

fsp_repeat_1_2_max_abs_diff=0.0
fsp_repeat_1_3_max_abs_diff=0.0
fsp_matrix_normalized_max_abs_diff=0.0
permuted_matrix_similarity_max_abs_diff=0.0

permuted_run.retcode=Success
permuted_run.raw_sum=-1.6787697354411202e-16
permuted_run.raw_residual_inf=1.159211204601055e-16
permuted_run.normalized_residual_inf=0.6905123318156889
permuted_vs_transpose_max_abs_diff=0.029096753221841617

Three repeated generated-RHS solves are byte-for-byte identical. A direct sparse-matrix solve produces exactly the same normalized result. The permuted sparse matrix satisfies the expected similarity transformation exactly, but its independent near-zero root points in a different roundoff direction, so post-hoc normalization produces a different invalid vector.

Normalization-constrained control

As a control, I replaced one redundant generator residual with the probability-mass equation before calling the same SSRootfind() algorithm:

function mass_constrained(f)
    return function (du, u, p, t)
        f(du, u, p, t)
        du[firstindex(du)] = sum(u) - 1
        return nothing
    end
end

No assertion or tolerance was changed. On the exact same parameter vector:

fsp_constrained_run_1.retcode=Success
fsp_constrained_run_1.raw_sum=0.9999999999999994
fsp_constrained_run_1.residual_max_abs=8.465450562766819e-16
fsp_constrained_run_1.marginal_a_max_abs_error=3.0531133177191805e-16
fsp_constrained_run_1.marginal_b_max_abs_error=2.498001805406602e-16
fsp_constrained_repeat_1_2_max_abs_diff=0.0

matrix_constrained_run.raw_sum=0.9999999999999994
matrix_constrained_run.residual_max_abs=8.396061623727746e-16
fsp_matrix_constrained_max_abs_diff=0.0

permuted_constrained_run.raw_sum=0.9999999999999998
permuted_constrained_run.residual_max_abs=7.325954079484553e-16
permuted_constrained_run.marginal_a_max_abs_error=3.191891195797325e-16
permuted_constrained_run.marginal_b_max_abs_error=3.0531133177191805e-16
permuted_constrained_vs_transpose_max_abs_diff=1.1107651257113993e-16

The reported residuals above are evaluated against the original, unconstrained generator RHS after the solve.

I also tried avoiding the root formulation with DynamicSS(Vern7()) at tight tolerances. It did not produce a result and was stopped after about 6m37s when its resident memory had grown to approximately 108 GiB, so that does not look like a practical test workaround.

History

• 594b482708fefd3464d6e49f2da41587b8265da1 (2022-03-10, “Added steady-state functionality”) introduced the homogeneous steady-state RHS/conversion and the test pattern solve(..., SSRootfind()); sol.u ./= sum(sol.u).
• 9ec1ebdaf8aed5f2460a49e140bf7921030bab55 (2024-07-12) added the species-permutation steady-state comparison and changed the truncation/random range.
• 90e969e25de24f181def691183985258c0cd36fe and b237d20a8b64d5dca093d38efb45d2fe01590732 later adjusted truncation and random parameter ranges for test stability, but did not make the homogeneous root problem well-posed.

Possible fix direction

The finite-state steady-state API needs an explicit normalization constraint, or a separate documented normalized-stationary problem/helper. Replacing one linearly dependent generator equation with sum(u) - 1 = 0 produced the expected stationary distribution for the generated RHS, sparse matrix, and permuted system without loosening any tests.

Because changing the public steady-state formulation may affect API semantics, I have not made a test-only workaround or silently changed the constructor here.

Start with the SteadyStateProblem and SSRootfind entry points, then inspect test/birthdeath2D.jl, especially the stationary assertions at lines 74-75 and 116. Review the historical steady-state formulation before choosing how normalization should be represented in the API. Done means a well-posed stationary solve that passes the existing marginal and permutation comparisons without weakened tolerances.