perf: avoid unnecessary exact elimination for triangular rational systems

[acgetchell]

Summary

Optimize exact rational solves for triangular and row-permuted triangular systems without losing the dense-system advantage of fraction-free elimination.

During the downstream adoption of RationalMatrix from #216 in acgetchell/delaunay#575, la-stack 0.4.6 was consistently slower than the removed rational Gaussian solver on a deterministic triangular family. The regression persists when matrix/vector construction is excluded. This is worth a focused optimization because the old implementation avoids forward-elimination arithmetic entirely on these systems, while the current backend performs cubic work. The observed end-to-end gap reaches roughly 11 µs per size-8 solve, and downstream exact LP certification can perform repeated basis solves.

These are synthetic triangular systems, not a measured distribution of real Delaunay LP bases or evidence that all sparse systems regress. No whole-triangulation speedup is claimed.

Measurements

Local aarch64 macOS; rustc 1.98.1 (48a229cea 2026-09-01); la-stack 0.4.6; Criterion 0.8.2. Both alternatives used a standalone perf profile inheriting release, with thin LTO, one codegen unit, and no debug information. Timed runs did not overlap compilation or other validation jobs.

Initial comparison, 20 samples with 1 second warm-up and 2 seconds measurement per case:

Triangular system size	Original Gaussian, including input copy	Full la-stack adapter	Change
2	0.738 µs	1.241 µs	+68%
4	2.981 µs	4.375 µs	+47%
6	6.703 µs	10.623 µs	+59%
8	12.266 µs	23.240 µs	+89%

A separate size-8 diagnostic run isolated construction, using 30 samples with 1 second warm-up and 3 seconds measurement:

Operation	Row-permuted triangular	Dense Hilbert control
Original Gaussian, including input copy	12.178 µs	114.290 µs
Full la-stack adapter	23.456 µs	35.123 µs
la-stack input construction only	3.109 µs	4.609 µs
RationalMatrix::solve with inputs constructed beforehand	20.858 µs	31.144 µs

These are Criterion central time estimates, not medians. Components were measured separately and include their own output destruction, so they need not sum exactly. Both solvers recovered the known exact solution before each benchmark. The baseline excludes Delaunay's removed rational-to-f64 round-trip probe, giving the legacy comparison an optimistic baseline.

Reproducer inputs and comparison

For each n = 2..=8, construct all entries with exact BigRational arithmetic and zero-based indices:

A[i,j] = (i + j + 1) / (3 + 2*i), if j >= i
         0,                       otherwise
rotate A's rows left by one
x_expected[j] = (j + 1) / 11
b = A * x_expected

The dense control uses A[i,j] = 1 / (i + j + 1) before applying the same row rotation and manufactured solution. Fractions must be formed as rationals, not floating-point division followed by conversion.

Use the captured original rational Gaussian solver as the reference. The legacy timed closure clones the matrix and RHS into its owned inputs. The upstream adapter uses try_with_rational_matrix!(n, ...), fills the zero matrix with set(row, col, value.clone()), builds the RHS with RationalVector::try_from_fn, calls solve, and collects the resulting array into a vector. Measure that adapter separately from solving a preconstructed RationalMatrix<8> and RationalVector<8>.

Suggested Criterion arguments after adding repository-owned benchmark cases:

--sample-size 30 --warm-up-time 1 --measurement-time 3

Diagnosis

After pivoting these triangular fixtures, every below-pivot coefficient is zero. The old Gaussian implementation skips the corresponding row update. It performs no forward-elimination arithmetic and proceeds to rational back-substitution.

The current RationalMatrix::solve clears row denominators and calls the shared integer backend. Its Bareiss loop updates every remaining row and column, including products with zero:

a[i,j] = (pivot * a[i,j] - a[i,k] * a[k,j]) / previous_pivot

For size 8 this executes 140 trailing-matrix updates plus 28 RHS updates. The zero elimination coefficient still leaves a pivot-ratio rescaling term. The public setters also re-canonicalize already-reduced inputs, accounting for some construction overhead, but the preconstructed-input experiment shows that construction is not the main gap. Timings do not separately apportion denominator clearing, forward elimination, and back-substitution.

Do not copy the old zero-coefficient continue directly into Bareiss: omitting the pivot-ratio scaling breaks its arithmetic invariant.

Proposed scope

• Investigate a structure-aware exact solve for triangular and row-permuted triangular matrices, with direct rational back-substitution over validated inputs. Retain the fraction-free path for general dense systems.
• Alternatively, use a mathematically justified sparse/fraction-free strategy that removes the measured overhead while preserving required scaling. Choose based on benchmark evidence rather than fixture-specific branches.
• Keep the optimization inside la-stack so downstream crates can continue removing duplicated generic elimination.
• Treat construction/canonicalization overhead as a separate concern. Do not bypass raw-denominator validation or expose unchecked rational storage to meet a timing target.
• Preserve the current public API and stable-Rust support; no new runtime dependency or unsafe code is needed for this scope.

Acceptance criteria

• Add permanent benchmarks for diagonal, upper/lower triangular, row-permuted triangular, general sparse, and dense rational systems through size 8, including both dyadic and non-dyadic coefficients. Separate construction from solve cost and keep fixture generation/correctness checks outside timed closures.
• On the reproduced triangular family, bring the full rational adapter to parity with or better than the original skip-zero Gaussian reference within measurement uncertainty. Report both kernel and end-to-end results, and retain the dense-system improvement without a material regression from structure detection.
• Check exact A * x == b and agreement with an independent reference across dimensions 0..=8, row permutations, mixed signs/denominators, zero RHS, and large numerators/denominators. Include cases that force the general fallback.
• Preserve typed singularity diagnostics for zero pivots, rank-deficient inputs, and duplicate rows, including cases encountered after row permutation. Preserve failure behavior for invalid raw rationals at public construction boundaries.
• Do not introduce implicit rounding, binary64 reconstruction, or weakened exactness/shape invariants. Keep arithmetic on the validated rational/integer domain.
• If shared Bareiss code changes, retain determinant value/sign, row-swap parity, and finite-f64 exact-solve regressions as well as rational-solve coverage.
• Run the repository's required validation and report before/after measurements on the same toolchain, profile, features, and fixtures.