Inline linear SCCs: 296×296 dense runtime solve for an ~18-unknown core; non-rank-tolerant; symbolic shrinkage intractable — solve sparse+rank-tolerant instead

[baggepinnen]

Summary

For multibody vehicle models, the inline-linear-SCC pass (DefaultReassembleAlgorithm(inline_linear_sccs = true)) emits a runtime dense A \ b over a system that is an order of magnitude larger than the mathematically-required core, and solves it with a non-rank-tolerant dense LU. Measured on a half-car suspension model (22 states after mtkcompile):

• The big algebraic SCC enters get_linear_scc_linsol with 396 rows.
• __reduce_linear_system! removes only 101 (the all-constant-coefficient rows) → a 296×296 dense runtime LU, evaluated every RHS call and during initialization reconstruction.
• Of the 396 rows, 384 are matched (torn) and 378 of those have a constant pivot coefficient — i.e. the irreducible simultaneous core is only ~18–24 unknowns (consistent with the same model compiled with inline_linear_sccs=false, which yields a 48-state DAE = 22 differential + 26 algebraic).
• The emitted 296×296 pattern is 1.6% dense and, after a Dulmage–Mendelsohn analysis, ~93% block-triangular given the tearing order; only the small core is truly coupled.

The inflation mechanism: rows whose non-pivot coefficients are symbolic (every product-rule-differentiated definitional equation — rotation-matrix entries, contact-frame vectors, etc.) fail the all(isconst, coeffs) gate in __reduce_linear_system! and are kept inside the numeric solve, even though tearing already matched them with constant pivots.

Why the obvious fix doesn't work (tested)

Relaxing the gate to "pivot coefficient constant" (which is sufficient for a division-free elimination, and the alias/substitution machinery in __reduce_linear_system!/get_new_mm handles symbolic alias coefficients as-is) shrinks the system but explodes symbolically:

1. First failure: iszero(::Num) inside get_new_mm's duplicate-entry pruning goes through fraction_iszero → expand, which is exponential for the deep product-of-sums coefficients → OutOfMemoryError. (Workaroundable with a cheap syntactic zero test — see suggestion below.)
2. With that bypassed, mtkcompile succeeds (~90 s) but ODEProblem construction stalls indefinitely: the substituted core coefficients are enormous expression trees and codegen blows up.

So full symbolic elimination of the torn rows is intractable; the inflation is an inherent artifact of numeric tearing. Sequential runtime evaluation of the eliminated rows is also impossible (they depend on the core unknowns — circular).

Suggested fixes

1. Solve the emitted system sparse and rank-tolerant at runtime. The pattern is already available (the ArrayMaker regions); a sparse LU/QR with fill-reducing ordering automatically exploits the ~93% triangular structure, so the 296-system costs about as much as the 18-core. Crucially it should be rank-revealing (or regularized): these acceleration/reaction systems are legitimately rank-deficient whenever constraint reactions are statically indeterminate (vehicle with ≥2 rigid suspension mounts on one body) — the current dense \ throws SingularException (LU) and the DyadCompilerPasses LDIV rewrite produces silent Inf (unpivoted QR). See JuliaComputing/DyadCompilerPasses.jl#40.
2. Independent robustness improvement: the duplicate-entry pruning in get_new_mm should use a cheap syntactic zero test (term construction already cancels syntactically identical terms; an entry that is mathematically zero but kept is harmless), e.g. an overridable prune_iszero(x) = iszero(x) specialized for Num to a syntactic check — exact iszero(::Num) can OOM via polynomial expansion.
3. Optionally expose the SCC composition (total/matched/const-pivot counts) as a debug aid — these numbers immediately reveal the inflation.

Context

Found while debugging SingularException(13) during initialization of suspension vehicle models (Multibody-style components). ModelingToolkit v11.26.8 (#master), ModelingToolkitTearing v1.14.1, StateSelection v1.9.3. Related: SciML/ModelingToolkit.jl#4237 (init reconstruction LU of indeterminate reactions), SciML/ModelingToolkit.jl#4607.

[baggepinnen]

Additional data point isolating what controls the emitted system size: compiling the same half-car with allow_symbolic = true (which grants solvable edges to the unknown-containing coefficients that otherwise fuse the SCC) reduces the state count (22→20; the quarter car goes 11→10 and still solves), but the emitted linear-system sizes are completely unchanged: [296, 15, 15, 13, 13, 12, 12] in both modes. So the inline system's size is governed by __reduce_linear_system!'s all-constant-coefficient gate, not by solvability semantics — consistent with the analysis above, and reinforcing that the practical fix is sparse + rank-tolerant runtime solving of the emitted block rather than changes to solvable-edge rules (whose symbolic elimination route we already showed explodes).

[oscardssmith]

Can you see if #97 helps at all?

[baggepinnen]

Data point from the HalfCar resolution (see #98 for the full identification): after the model-level state-priority fixes (JuliaComputing/MultibodyComponents.jl#67), all small inline-linear blocks are full rank at the real parameter point and the model integrates end-to-end, so on this model #95 is now purely a performance concern (the ~300×300 dense LU per RHS call; it was already full rank per the measurements above), no longer entangled with rank-deficiency. The runtime probe used for these measurements is committed as test/probe_halfcar_blocks.jl on the MultibodyComponents suspension branch.

[baggepinnen]

More evidence for the pivot-criterion + sequential-emission direction, from the half/full car models (companion measurement to the Robot6DOF investigation; baseline = MultibodyComponents suspension branch, fbc/halfcar-stack StateSelection stack).

First, the model-side check is clean here: unlike Robot6DOF (isroot=false), the cars already use the MSL sweep direction everywhere (MBC Revolute default rooted=FrameA = kinematics outward, forces inward). I swept all 8 symmetric rooted flips over the cars' revolutes (heave/roll chain, four-bar r2, wheel spin): the big SCC is insensitive (295–298 emitted in every config) and every flip regresses — flipping the roll revolute spawns seven 3×3 position-level inverse-rotation solves, and flipping r2/wheel_rotation explodes the state count (14 → 34/52, rotation-matrix entries become states). So the numbers below are from a correctly-directed model; the "provided the model is correctly implemented" caveat is satisfied.

Pivot-criterion headroom, measured in __reduce_linear_system! (ENV-gated count of matched / all-const-coefficient / const-pivot rows + torn vars, inserted just before the elimination loop):

	SCC size N	matched in SCC	all-const rows (current gate)	const-pivot rows	torn core
HalfCar big SCC	396	384	99 → emits 297	379 → would emit ~17	12
FullCar big SCC	791	741	223 → emits 571	729 → would emit ~62	50
per-suspension SCC (×2/×4)	47	44	1 → emits 46	44 → would emit 3	3

LU-flop ratios: ≈5000× (HalfCar, 297³→17³) and ≈780× (FullCar, 571³→62³). For scale: FullCar currently integrates at ~0.7 s per accepted step with Rodas5P, essentially all of it in the emitted dense LU + FD/AD Jacobians thereof, so this is where the entire solve cost lives.

Re "shrinks the system but explodes symbolically" (top post): that was full symbolic elimination of relaxed-gate rows. With the residual numeric systems this small (~17/~62 — only the genuinely non-const-pivot rows), the Robot6DOF note's alternative — emit const-pivot matched rows as sequential assignments and solve only the small core numerically (the OM strategy) — avoids the fill-in entirely and would land these models near their torn cores (12/50). The remaining requirement is a nonzero guard for parameter pivots (maybe_zeros-style), which on these models is what distinguishes the 379 const-pivot rows from the 384 matched ones.

Measurement scripts/results: rooted_sweep_results.txt, scc_stats_results.txt on the MultibodyComponents checkout; the stats instrumentation is ~20 lines in __reduce_linear_system!, happy to push it ENV-gated if useful.

[baggepinnen]

Summary

The emitted inline-linear-SCC runtime solve A \ b is not rank-tolerant, and on multibody vehicle models the acceleration block is numerically singular / severely ill-conditioned at the model's axis-aligned joint configurations. There is currently no runtime-solve option that is both correct and robust there:

runtime solve	behavior on the large HalfCar block at a consistent state
DyadCompilerPasses LDIV_RULE (optimize=:basic, unpivoted)	silently wrong, finite — a mirrored suspension wheel's spin acceleration comes out 0.117 instead of ~0
plain pivoted \ (optimize=:none)	throws ArgumentError: matrix contains Infs or NaNs on the exactly-singular step during integration
rank-tolerant (pivoted-QR / pinv min-norm)	correct (~0)

So optimize=:basic integrates to a wrong trajectory and optimize=:none crashes; only a rank-tolerant solve is both correct and non-throwing. This is the runtime-solve half of #95 (the "non-rank-tolerant" point), separated out with a concrete repro and a clear contrast against the construction side.

Evidence the block, not the reduction, is at fault

• The inline-linear reduction is exact: reconstructing the full 396-row HalfCar acceleration block from the emitted reduced (A_red, b_red, aliases, constants) reproduces the ground-truth solution (full-block residual ~1e-10). Pivot-solvability inline-linear row elimination with size-budgeted composition (#95) #105 confirms the same for the pivot-gate reduction.
• Solving the unreduced 396-row block directly with a pivoted/pinv solve at the consistent state gives the mirrored-wheel value ~0 (the block is full-rank there but very ill-conditioned; at later integration steps it becomes exactly singular — hence the :none throw).
• The wrong 0.117 is therefore produced purely by the unpivoted LDIV solve, not by any reduction.

Where

INLINE_LINEAR_SCC_OP = (\) (lib/ModelingToolkitTearing/src/reassemble.jl) and the DyadCompilerPasses LDIV_RULE rewrite of that \ (get_factorization uses unpivoted qr(A); JuliaComputing/DyadCompilerPasses.jl#40/#41). A column-pivoted/rank-revealing factorization (qr(A, ColumnNorm())) or a min-norm pseudo-solve would make the large block correct; the legitimately rank-deficient consistent blocks (the static-reaction-indeterminate ones) then get a valid min-norm solution instead of garbage or a throw.

Mitigation available today

#105 shrinks the emitted block (HalfCar 297→64) enough that it is no longer ill-conditioned for these models, so all three solve variants agree and the cars integrate correctly — but that is model-dependent and does not fix the underlying non-rank-tolerant runtime solve; a differently-conditioned block elsewhere would still hit this.

Repro

MultibodyComponents.examples.suspension.HalfCar on the #99+#100+#103 stack; build multibody(hc), init at the canonical defs, then evaluate the RHS with optimize=:basic (→ 0.117 for excited_suspension_l.wheel_rotation.w's derivative) vs optimize=:none (→ correct at the init point, throws mid-solve). Probe scripts available on request.

[AayushSabharwal]

JuliaSymbolics/Symbolics.jl#1893

Should address iszero being slow.