Add analyze_initialization_jacobian — diagnose rank-deficient initialization

[baggepinnen]

Motivation

Rank-deficient initialization systems are a recurring source of intermittent ("flaky") initialization failures. When the init Jacobian doesn't have full rank, the solve can land on different solutions — or fail — depending on the solver and on the (nondeterministic) ordering of the assembled init unknowns/equations: some orderings hit a singular realization and throw SingularException, others succeed, so the same model fails maybe 1 run in 10.

MTK already warns structurally (overdetermined/underdetermined … N equations for M unknowns, plus singular_check's problematic-variable list), but it's often hard to tell which unknowns form the free direction, and — when an init has more equations than unknowns yet is still underdetermined — which equations are redundant. I hit this in a multibody model and ended up computing the init-Jacobian SVD by hand; it seems generally useful.

What this adds

analyze_initialization_jacobian(prob; rtol, atol, threshold, verbose) evaluates the Jacobian of the initialization residual at the initial guess u0 (via ForwardDiff), computes its SVD, and reports both sides of the rank deficiency:

• Underdetermined unknowns — the right null space (directions the init is free to move along). Pinning these makes the init well-posed and deterministic.
• Redundant equations — the left null space (combinations of init equations whose Jacobian rows are linearly dependent at u0, so they don't locally constrain any further DOF). This explains the confusing case of "more equations than unknowns, yet still underdetermined."

prob may carry initialization data (an ODEProblem/DAEProblem built from a System) or be an initialization NonlinearProblem/NonlinearLeastSquaresProblem directly.

Per-item weights are the diagonal of the respective null-space projector (squared row norm over an orthonormal basis), so they're invariant to the arbitrary basis chosen within the null space. Returns (; jacobian, singular_values, rank, nullity, redundancy, underdetermined_unknowns, redundant_equations) and prints a report by default.

Example

julia> analyze_initialization_jacobian(prob)
Initialization Jacobian rank analysis
  residual Jacobian: 58×51, rank ≈ 47, nullity ≈ 4, redundancy ≈ 11
  smallest 5 singular value(s): [0.41, 0.2, 0.07, 3e-15, 1e-15]
  Underdetermined unknowns (participation ∈ [0, 1]):
      0.746  rotational_losses_bl₊w_rel(t)
      ...
  Redundant equations (participation ∈ [0, 1]):
      0.62   0 ~ √(…) - springdamper₊length(t)
      ...

which immediately points at both the unpinned degrees of freedom to constrain and the equations that aren't pulling their weight.

Notes

• Added to src/initialization.jl, exported, with a docstring and tests in test/initialization_jacobian_analysis.jl (registered under the Initialization group): an underdetermined DAE (nonzero nullity), an overdetermined system with a redundant init equation (nonzero redundancy, equation named), a determined algebraic init (full rank), and a fully-IC-determined system (no init problem — handled gracefully).
• The Jacobian is evaluated at a single point, so it reports local rank structure (the docstring calls this out).
• Happy to rename or fold into existing init diagnostics if you'd prefer.
• I couldn't run the full MTK test suite in my (offline) environment; the function's logic and all four test cases were validated against MTK on small models + a downstream multibody model. CI here will exercise the added testset.

🤖 Generated with Claude Code