Mixed Precision on MTL5

Status: in progress

The equivalence gate for the port is in place. All five families run on MTL5. Four are bit-identical: the products and sums, the Cholesky factor with its SPD solves, the Householder QR with the MSCKF null-space projection, and the square-root form’s triangular solves. The fifth, the SfM eigen/SVD family, runs on MTL5 5.20’s eigen_symmetric_jacobi. It cannot keep its bits, but it gives the same decisions and a more accurate null vector. The issues: #448 (MTL5 port), #449 (genericity over non-IEEE types), #450 (per-stage arithmetic and precision characterization) and #451 (the MTL5 gap register). The benches and the type-scaled invariants they rely on are in place, and every stage bench measures the bits it loses in each type (--precision, see below). The MSCKF backend can run each stage in its own type (per-stage types).

The same effort that makes every stage observable also puts the pipeline’s linear algebra on MTL5, generic in the arithmetic type. Running the pipeline in several number systems then becomes a diagnostic of its own: a stage whose output changes materially between double, posit<32,2> and a wider type is numerically fragile.

Why it matters here

Several of the pipeline’s failures were numerical:

• the square-root covariance diverging on real data;
• a non-finite δx from a barely-conditioned update;
• the Clang-only divergence of #437.

With per-stage arithmetic, “is this algorithmic or numerical?” becomes a measurement. Run the stage in a wider and a narrower type on the same fixture, and see whether its invariants move.

Today vs the port

Before the port the pipeline made no mtl:: calls. Its only MTL5 touchpoint was the mtl::Scalar concept behind math::Scalar, and its linear algebra was hand-rolled. DynMat stays as the storage, and each family below moves its operations onto MTL5 views of it:

Hand-rolled	MTL5 replacement	Status
msckf::DynMat<T>, a row-major std::vector holding P, F, G, H, S, K	zero-copy mtl::mat::dense2D<T> views	done
GEMM, add	operation/mult.hpp (dense, and CRS × dense for sparse left operands), mat/operators.hpp	done
Cholesky, SPD solve	cholesky_factor (term-by-term order), two trsm calls	done
GEMV, transpose	operation/mult.hpp, trans.hpp	the square-root gain product Kbar·y is mult; transpose stays a copy
A·P·Aᵀ sandwiches (clone, predict, S, Joseph)	mult (syrk/symm don’t fit, see below)	on mult through mul; a congruence primitive is a gap
Householder QR (R only)	qr_factor (textbook reflector)	done
Triangular solves (square-root form)	trsm over a map_view block	done
MSCKF left null space	householder + apply_householder_left (textbook reflector)	done
Symmetric eigen, 3×3 SVD	eigen_symmetric_jacobi (MTL5 5.20); svd3 builds on it	done, not bit-identical (see below)
Principal-submatrix gather	mat/view/map_view.hpp	views where an MTL5 kernel reads the block; plain copies stay loops

How the port is gated

1. Safety gate first. At T = double, the MTL5 port must reproduce today’s results: the same EuRoC trajectories to ≤ 1e-9 relative, and all current tests, gates and inspector outputs unchanged. The bit-identical equivalence harness used for the stage refactor is the model. One family is exempt from bit identity: the SfM eigensolver, whose MTL5 kernel rounds differently. Its six sfm_* bench outputs are gated on the same decisions and on null-vector accuracy instead (see the eigen family).

2. Stage by stage, each behind its bench. A stage is ported only when its bench passes on its fixtures in every type. A whole-pipeline check closes the port.

3. Fix the sites that break or narrow a non-IEEE T:

   • non-ADL math (std::sqrt, std::acos on T);
   • static_cast<double> on the estimation path;
   • fixed ridges and epsilons like T(1e-9) that vanish in 16-bit types, which should be T-relative;
   • the double timestamp mixed into T.

   By design, timestamps stay double and the telemetry accumulators in eval/ stay double.

4. Per-stage arithmetic. Covariance and update math can run in one type while Jacobian assembly runs in another. A precision matrix of stages × types becomes the characterization report, and the diagnosis of #437.

5. No latency regression. vio_bench per-frame latency (p50, p99) on EuRoC V1_01 must stay at its pre-port level, or the change must explain the difference.

Genericity, checked by a strict type (#449)

Point 3’s sites are found by the compiler. StrictDouble (tools/bench/strict_double.hpp) is an IEEE double whose conversions out (to double, float or an integer) are all explicit. Construction from an arithmetic type stays implicit, as for a posit. Every stage bench instantiates its stage in it, so a double x = t;, a std::sqrt(t) that picks the double overload, or a T passed where a double is taken, fails to compile. Every operation is the double operation, so a generic stage computes exactly what it computes in double. A test holds every bench to that:

./build/tools/s6e_ekf_update_bench --types double,strict_double --variant all
#   max |output - double output| = 0     (every fixture, every variant)

The first run compiled cleanly and matched double bit for bit on every fixture and variant of all 14 benches, the C_filter_loop composition included, which runs the whole estimator. The non-ADL math and static_cast<double> sites the audit listed had already gone: Universal’s posits convert to double only explicitly too, so the posit32 benches removed them. What strict_double adds is a check at full double precision that genericity costs nothing at T = double, and a compile error the moment a narrowing comes back. The static_cast<double> that remain on the path (S1’s diagnostics, the NIS telemetry, the double timestamp) are explicit and are telemetry or timestamps by design.

The type is left out on MSVC, whose Release optimizer crashes (C1001) compiling the stages with it. The check runs on gcc and clang, and the posit types already enforce explicit conversions on MSVC.

Thresholds that scale with T, and what stays IEEE (#449)

The audit listed the fixed T(1e-9) ridges and SO(3)‘s 1e-5 small-angle threshold as constants that vanish in a coarse type. Measured in the benches, they split into two kinds:

• Tolerances scale with T. The small-angle hand-off in SO(3) and Sim(3) is lie::detail::small_angle<T>() = max(1e-5, ε_T^(1/3)). In float, 1 − cos θ is already zero at θ = 1e-4, so the closed-form Jacobian lost its W/2 term there. Below ε^(1/3) the series’ truncation error is under ε, and above it the closed forms lose at most ε^(2/3). Now the float Jacobians at θ = 1e-4 match double to 1e-6. Sim(3)‘s scale-rate integrals $\int_0^1 u^k e^{\sigma u}\,du$ cancel like σᵏ, so no shared hand-off serves them. They are now a series summed to convergence for |σ| < 1, accurate to about an ulp in every type. Sim(3) is not on the VIO path. Two more tolerances use math::at_least_roundoff(nominal, scale) = max(nominal, 16·ε_T·|scale|):

  • the PSD check’s ridge, floored at 16 ulps of ‖P‖ (its largest diagonal entry), so that rounding at float’s level is not read as a negative eigenvalue. A floor per diagonal entry was not enough: the clone copies make P exactly singular, and rounding puts that zero eigenvalue near −ε·‖P‖, which on a low-variance row is far below ε times its own diagonal. An end-to-end float run hits this within seven frames. The check is a diagnostic (probes and inspectors); no estimation decision reads it;
  • the initializer’s up-vector degeneracy test, which compares a sine in [0, 1].

• Regularizing ridges stay fixed. Triangulation’s Gauss-Newton ridge and the initializer’s gyro-bias and velocity ridges guard a zero diagonal, an unexcited unknown. Every bench type represents 1e-9 there: double, float and posit32 hold it, and posit16 rounds it up to its minpos (3.7e-9), not to zero. IEEE half and Q16.16 fixed point round it to 0; a port to those needs a rescaled problem. A T-scaled ridge was tried and measured as harmful:

  • floored at 16·ε_T, the velocity/gravity/scale alignment’s float scale error in the S1 SfM bench went from rounding level to 17.2 %;
  • floored at 16 ulps of each diagonal entry, it went to 24.8 %.

  A ridge λ shrinks the solution along each eigenvector of the normal matrix by λ/(σᵢ + λ). The alignment’s scale-carrying mode has σ = 5.5e-6, so 1.9e-6 removes about a quarter of it. These ridges stay at 1e-9 in every type.

For double (and long double) small_angle returns 1e-5 exactly. at_least_roundoff returns the old 1e-9 unless the entry it scales with exceeds about 2.8e5, where 16·ε·|scale| passes 1e-9; that is a property of the input, not a guarantee. As measured, the bench snapshot and the EuRoC trajectories are unchanged bit for bit. The error models, the measurements behind each choice and guidance for hardware formats are on Numerical Thresholds Across Arithmetic Types.

What stays IEEE by design, outside T:

Part	Type	Why
KLT tracker (cv/klt.hpp)	double internally, float points	an image-domain front end; pixels are the measurement, not the estimate
VioConfig	double	a parsed configuration POD (see the configuration discipline); cast to T once when the backend initializes
Timestamps	double	time is an index, not a state; dt is cast to T for the propagation step, and the backends then reset the clock to the double sample time, so a narrow T’s rounded Σ T(dt) never reaches the reported pose times
eval/ telemetry (NEES, NIS, innovation statistics)	double	diagnostics accumulate over a whole run; they read the estimate, never feed it

End to end in every type (#449)

The C_filter_loop composition runs the whole backend loop (S2 propagation through S9 marginalization) in each type. The estimator itself, front end and backend together, also runs end to end on a synthetic scene whose IMU and images agree: a level body accelerating at 2 m/s² under a plane of landmarks 3 m away, projected through the estimator’s pinhole camera, so triangulation sees parallax and every S6 update passes the gate. The front end stays IEEE inside (see the table above) and hands T-valued observations to a backend that runs entirely in T.

• tests/sdk/vio_estimator.cpp runs it in float, posit<32,2> and long double. It requires accepted updates, and positions within an absolute 1000·max(ε_T, ε_double) m of the double run (1000 ulps at a 1 m scale; double’s ε bounds long double, whose deviation is double’s own rounding). Measured over 16 cm: float 1.4e-5 m (bound 1.2e-4), posit32 5.4e-7 m (bound 7.5e-6), long double 9.0e-15 m (bound 2.2e-13).
• tests/tools/strict_double.cpp runs it in strict_double and requires every pose to match double’s bits.

Every invariant passes in every type. Worst |output − double output| on the loop’s fixtures:

Fixture	float	posit32	long double
known_answer_constant_velocity	1.0e-7	4.7e-9	1.1e-16
ground_truth_synthetic_three_frames	7.9e-5	1.1e-5	1.4e-13
captured_synthetic_three_frames	9.9e-4	7.6e-5	1.7e-12

posit32 lands 7–22× closer to double than float on the loop. It also runs about 1,300× slower than float on the three-frame fixtures, because it is emulated in software. S4 (the KLT front end) is the one bench that runs in double only, by the design above.

The equivalence gate

Two tools hold a port to “same results”. Run each before the port to take a reference, and after it to compare:

# Stage benches: every built-in fixture, every variant, in double (132 outputs).
./build/tools/bench_snapshot --write before.json              # before the port
./build/tools/bench_snapshot --compare before.json            # after: 1e-12 normwise, per output

# EuRoC: every tracked pose and its σ, at full precision, per sequence.
./build/tools/vio_trajectory_report --dataset …/V1_01_easy/mav0 --label V1_01 --out ref/    # before
./build/tools/vio_trajectory_report --dataset …/V1_01_easy/mav0 --label V1_01 --out after/  # after
./build/tools/bench_snapshot --compare-trajectory ref/V1_01_trajectory.csv after/V1_01_trajectory.csv
                                                               # the gate: 1e-9 normwise, per column
cmp ref/V1_01_trajectory.csv after/V1_01_trajectory.csv       # optional: bit-identical too?

Both compares fail on anything they can’t check:

• a fixture that throws;
• an output or column present on one side only;
• a non-finite value that isn’t the same kind on both sides (NaN, or an infinity of the same sign).

What the first family showed

The products port reproduces every bench output and all three EuRoC trajectories (V1_01, MH_05, V2_03) bit for bit: MTL5 builds each entry in the same k order as the hand-rolled loop. It also cost +7 % p50 and +17 % p99 per-frame latency.

• MTL5 5.12’s generic dense product is an i-j-k dot product. It walks B down a column and can’t vectorize, so it ran 1.5–2.6× slower than the i-k-j loop it replaced at the filter’s sizes (81×81, 81×15).
• MTL5 has no CRS × dense-matrix product. The hand-rolled loop skipped zero terms, which made F·P (F is the identity on the clone rows) and (I − KH)·P cheap.

Both were fixed upstream rather than worked around: mtl5#549 is a row-major i-k-j kernel, bit-identical by construction, and mtl5#550 is mult(compressed2D, dense2D, dense2D), both in mtl5#552, released in MTL5 5.14. With it the port is within 1 % of baseline p50 latency (+0.7 % p50 and +1.4 % p99 for both families together: 30.79 → 31.02 ms and 37.32 → 37.84 ms, medians of three interleaved runs on V1_01), with byte-identical trajectories. What it still pays for is building a CRS copy per product, which a non-owning CRS view in MTL5 would remove.

What the second family showed

The second family is the Cholesky factor and the SPD solves behind it. They serve the update’s S⁻¹ and Kalman gain, the χ² gate, initialization, triangulation and the NEES checks.

• The solves port bit for bit. cholesky_solve is now two mtl::trsm calls, forward through L and back through trans(L). trsm’s generic path subtracts each term in the same order as the hand-rolled substitution. All 117 bench outputs and all three EuRoC trajectories are identical, and latency is unchanged (V1_01 p50 30.86 → 30.85 ms, p99 37.51 → 37.49 ms, three interleaved runs each).
• The factor does not. MTL5’s cholesky_factor forms each entry as a − (Σ products). The hand-rolled factor subtracts term by term. Both are backward stable, but they round differently: 37 of 117 bench outputs move, and the filter loop moves by up to 2.5e-11.
• On EuRoC, that rounding change is a different run. On V2_03, 2 of 16,265 update accept/reject decisions flip, and the trajectory goes a different way from there (ATE 0.269 → 0.261 m, worst pose about 15 cm apart). The filter’s discrete decisions amplify any rounding change. So the ≤ 1e-9 trajectory criterion only holds for a port that keeps the bits; a port that rounds differently needs a statistical comparison instead.
• Neither ordering is more accurate. In backward error ‖A − LLᵀ‖/‖A‖ at n = 81, term by term is about 2× better on ill-conditioned spectra like the filter covariance (cond 1e6–1e10), and the MTL5 order is about 1.4× better on well-conditioned ones.

So the factor waited for MTL5. mtl5#556 added an opt-in order, cholesky_factor(A, mtl::tag::subtract_term_by_term{}), which starts from the matrix entry and subtracts each product in turn. The same change made the pivot test reject a NaN, as cortex’s factor did. With that order the factor ports bit for bit as well.

QR and the null space, measured before porting

MTL5’s householder is the LAPACK reflector: scaled, v(0) = 1, with the sign chosen so R’s diagonal is positive. cortex’s QR and MSCKF null-space projection use the textbook reflector (unscaled v = x − αe₁, α = −sign(x₀)‖x‖). A prototype port showed:

• V1_01 and MH_05 summary metrics identical to 4 digits, but the trajectories at 1.2e-8 (outside the 1e-9 gate);
• one of 16,265 update decisions flipping on V2_03, after which the run diverges;
• a different, equally valid, null-space basis from S6b.

As with the Cholesky factor, MTL5 was asked for an opt-in textbook reflector (mtl5#559). MTL5 5.15 has it as mtl::tag::textbook_reflector, taken by householder, apply_householder_left and qr_factor:

// R of A (qr.hpp): MTL5 overwrites V with R on and above the diagonal
mtl::qr_factor(V, tau, mtl::tag::textbook_reflector{});

// MSCKF null-space projection (msckf_nullspace.hpp), column c of H_f
auto h = mtl::householder(x, mtl::tag::textbook_reflector{});
mtl::apply_householder_left(M, h.v, h.vtv, c, c + 1, mtl::tag::textbook_reflector{});

With the tag, the port is bit for bit: all 126 bench outputs at relative 0, and the V1_01, MH_05 and V2_03 trajectories byte-identical. A column whose norm is NaN is now skipped in the projection, as it already was in the QR, instead of being reflected.

The square-root form’s solves

SqrtCovariance forward-solves the innovation factor Re twice: in the χ² gate (γ = ‖Re⁻¹r‖²) and in the array update (δx = Kbar·Re⁻¹r). Both now call mtl::trsm on a map_view of Re’s block inside the triangularized array, so nothing is copied. trsm was chosen over lower_trisolve and trsv because of the order in which it sums:

• trsm’s generic path subtracts each product from b in turn and divides by the pivot. That is the hand-rolled loop’s order, so the results are bit-identical.
• lower_trisolve and trsv subtract one accumulated sum, b − Σ. That rounds differently, the same difference as the Cholesky orders above.

A zero pivot still stops the solve where it did. The gate rejects but reports γ over the rows solved before the pivot, which the NIS sample records. The update is a no-op. The solve never divides by the zero, which matters for types whose division by zero throws. Kbar·y is mtl::mult over another block view, and its generic product sums in ascending order from zero, as the loop did.

The bench gate covers the array update (the S6e sqrt_array variant), but the snapshot never calls the square-root gate. So the port was also checked against the old header directly: 4,000 random systems, 537 of them rejected at a zero pivot, with γ, δx and S identical bit for bit on gcc and clang.

The eigen and SVD family: a convergence bug first

The symmetric eigensolver is used only in the SfM front end of dynamic initialization: the 8-point null space (N = 9), triangulation (N = 4), the DLT PnP (N = 12) and svd3. No MTL5 routine reproduces cortex’s cyclic Jacobi, so this family cannot keep its bits. EuRoC trajectories stay byte-identical whatever it does, because dynamic init is attempted there but never succeeds. Measuring the candidates turned up a bug: cortex’s Jacobi stopped only at an exactly zero off-diagonal, which rounding never produces, so every call ran all 100 sweeps (#497). Null-vector accuracy and time per call on AᵀA systems with a known null vector (N = 9):

Solver	Null-vector angle	Time per call
cortex Jacobi, 100 sweeps (before #497)	3.3e-9	164 µs
cortex Jacobi, converged (off ≤ ε²‖A‖²_F)	3.2e-9	7 µs
mtl::eigen_symmetric (tridiagonal QR), tol = ε	9.0e-9	4 µs
mtl::svd of A (one-sided Jacobi)	3.2e-9	69 µs

Fixing the stopping rule made the solver 22× faster with the same accuracy. One window build in the backend dropped from about 0.55 s to about 28 ms (Release).

The SfM benches of #472 were in place by then, and none of MTL5’s kernels at the time was a drop-in for this family. eigen_symmetric was fast but less accurate on the null vector, and svd of A was accurate but 10× slower than the converged Jacobi. The missing kernels were filed upstream under mtl5#567:

• a cyclic Jacobi symmetric eigensolver (mtl5#564);
• deflation and an ε-scaled default tolerance in eigen_symmetric (mtl5#565);
• a V-only mode and QR preconditioning for svd (mtl5#566).

All three shipped, and MTL5 5.20 is pinned.

The eigen family on eigen_symmetric_jacobi

sfm::jacobi_eigh is now a thin wrapper over mtl::eigen_symmetric_jacobi. Its eigenvalues come back ascending. Every caller takes the argmin, and svd3 sorts on its own, so no call site changed. svd3 stays a 3×3 built from the eigendecomposition of MᵀM, now on the MTL5 solver: svd of A would cost 10× more for no accuracy the callers use.

MTL5’s Jacobi applies the same rotation as cortex’s did, but rounds differently in two places:

• Stopping rule. It skips a pair once |a_pq| ≤ ε·√|a_pp|·√|a_qq| (Demmel & Veselić) and stops after a sweep that rotates nothing. cortex stopped once the off-diagonal mass was below ε²‖A‖²_F.
• Update. It updates through τ = s/(1+c) instead of c and s directly.

So this family cannot keep its bits. It is gated on what it decides and on how accurate it is.

The relative stopping rule is a real accuracy gain. Null-vector error ‖v ∓ z‖ against an exact null vector z, on AᵀA with A = 3N×N, rows orthogonal to a known vector, then columns scaled over the given range. 500 trials per row, worst case (Release):

Column scale	N = 4: old / MTL5	N = 9: old / MTL5	N = 12: old / MTL5
1	1.1e-15 / 1.2e-15	1.9e-15 / 1.3e-15	2.7e-15 / 1.9e-15
10⁴	1.5e-9 / 1.7e-13	2.2e-8 / 3.1e-14	2.4e-8 / 1.8e-13
10⁶	3.9e-8 / 3.9e-12	2.9e-4 / 2.4e-13	2.4e-4 / 1.4e-12
10⁸	3.4e-1 / 1.4e-12	3.2e-1 / 1.1e-12	7.1e-1 / 2.8e-12

The Frobenius stop treats an entry below ε‖A‖ as converged even when it still couples the two smallest eigenvalues. Once the scales spread, that loses the null vector. MTL5’s stop does not. On well-scaled matrices both are at the rounding floor. A call costs 5–20% more (about 7 µs → 7.3 µs at N = 9).

The gate:

• Bench snapshot: 126 of 132 outputs are bit-identical. The 6 that change are the sfm_* variants of S1, the only callers, by at most 1.8e-12 relative (normwise).
• Captured EuRoC SfM windows:
  • MH_05 at rest (PnP fails) and moving with two-view failing replay bit for bit.
  • MH_05 aligned keeps its outcome (window built, alignment succeeds). Its outputs move by at most 1.7e-11 absolute, which is rounding-level change in the null vectors carried through RANSAC, PnP and the visual-inertial alignment. Its expected output was re-recorded from its unchanged input.

The 1e-12 compare above rejects those six outputs, so the exception has its own procedure, run against a snapshot taken before the change:

# 1. Only the exempt outputs moved: prints nothing.
./build/tools/bench_snapshot --compare before.json --rel 0 | grep DIFFERS \
    | grep -v '^ *DIFFERS  S1_initialization/[^/ ]*/sfm_'
# 2. ...and only by rounding: passes (worst 1.8e-12).
./build/tools/bench_snapshot --compare before.json --rel 1e-11
# 3. Same decisions, and the null vector kept on a badly scaled matrix.
ctest --test-dir build -R "S1 SfM|jacobi_eigh"

Step 3 covers the decisions: the captured EuRoC windows stop where they did, and the aligned one replays its re-recorded output bit for bit. It also covers the accuracy test. The exception covers this change only. A snapshot taken after it is the reference for the next port, and the plain 1e-12 compare applies again.

Why syrk and symm don’t carry the sandwiches

MTL5 has syrk, symm and syr2k, but none of them is the A·P·Aᵀ congruence the covariance needs:

• syrk forms A·Aᵀ. Using it means factoring P = LLᵀ first, which is the square-root form, already on MTL5.
• symm forms P·B for a symmetric P. Its generic path is a full GEMM that reads both triangles, and its BLAS path is for column-major operands only, while DynMat’s views are row-major. It saves nothing over mult, and it would compute (G·P)·Gᵀ as G·(P·Gᵀ), which rounds differently.
• F and G are sparse. mul already takes them through the CRS × dense product, which beats either dense kernel.

What would help is a triangular-output product, BLAS gemmt or a fused A·P·Aᵀ that computes one triangle and mirrors it. That halves the second product and makes symmetrize unnecessary. It is filed as mtl5#656. Adopting it won’t be bit-identical: averaging the two triangles rounds differently from mirroring one.

Per-stage precision characterization (#450)

Every stage bench has a --precision mode. It runs the stage once per arithmetic type (double, float, posit<32,2>, posit<16,1> and long double) on the bench’s fixtures and on a conditioning ladder, a short series of inputs along the axis that makes the stage hard: the condition number of P, parallax, depth, measurement σ, step count. For each output group (attitudes, positions, P, δx, NIS …) it reports how many of the type’s bits the stage lost:

./build/tools/s6e_ekf_update_bench --precision --csv out/   # → out/S6e_ekf_update_precision.csv

  ladder_sigma_0.05px [known_answer]  variant=shipped  log10_cond_P=15.65  log10_cond_S=4.279
    float        dx 0.0027 14.5/14.7       P_plus 0.00024 11.0/11.5  invariants 10/10
    posit32      dx 0.00028 15.2/15.3      P_plus 1.6e-05 11.1/11.4  invariants 10/10
    posit16      dx 1 12.0/0.0             P_plus 0.96 11.9/0.0      invariants 1/2  first FAIL innovation.spd
    long_double  dx 3e-12 13.7/13.7        P_plus 1e-12 12.2/12.2    invariants 10/10

Each cell is <group> <relative error> <bits lost>/<stage bits lost>:

• Relative error. The error is ‖y_T − y_double‖∞ / max(‖y_double‖∞, floor). The floor is the group’s noise level: half a pixel for a residual, 1 for an NIS, 1 µ for δx. Without it, a quantity that is about zero in double, such as a residual at the true poses, would count as having lost every bit.
• Bits lost. This is log₂(relative error / ε_T), capped at T’s precision. It is the stage’s error amplification at that input, an empirical condition number.
• Stage bits lost. This measures the same thing against double run on the input as T holds it. It separates the stage’s own arithmetic from the input’s rounding into T. Total and stage agree when the stage, not the input, is the source.
• Long double. The reference is the double run, so a long double row measures double’s own error, counted in ε_double.

The invariants run in every type, so a type that breaks the stage’s contract shows up as a failure. The mode is a characterization: it exits non-zero only if the double reference itself fails.

What the stages lose

The worst bits lost over each stage’s fixtures and ladder, of the 23 a float carries, 27 for posit32 near 1, 12 for posit16, and 52 for double. “fail” counts the fixtures on which an invariant broke. “flip” counts the fixtures on which a decision changed: a feature resolved or not, an update accepted or not.

Stage	float	posit32	posit16	long double (double’s own)	What drives it
S0 sensor model	1.4	2.7	9.8, 1/9 fail	0.5	distortion at the image corner; posit32’s loss is the pixel coordinates’ rounding (hundreds, where its taper leaves fewer bits); posit16’s finite-difference Jacobian fails at the corner
S1 initialization	2.4	3.4	6.0	2.0	the gravity and bias estimates over the window
S2 propagation	5.9	6.6	9.0	5.4	step count: 200 steps cost about 5 bits in every type, 25 steps about 2
S3 augmentation	0.0	0.0	0.8	0.0	none: block copies, exact at κ(P) from 10⁰ to 10¹²
S5 triangulation	2.9	3.7	12.0, 2/7 fail	3.8	pixel noise at low parallax (noise-free, even 0.25° costs float nothing); posit16 cannot resolve a feature at 1° or less
S6a Jacobians	9.1 (r)	8.8 (r)	9.8 (r)	9.0 (r)	the residual z − h(x) cancels: its error is ε·|z| in every type, 9 bits of the half-pixel floor
S6b null space	4.1	5.0	8.1	4.0	depth; posit32 loses more at 64 m, where H_f’s entries are far from 1
S6c compression	0.0	0.0	1.8	0.0	none: QR is backward stable
S6d gating	3.4	6.7	12.0, 2/6 fail	6.0	κ(S), dof; posit16’s S is not SPD at 15 rows and up
S6e EKF update	14.5	15.3	12.0, 3/7 fail	13.7	κ(S): δx loses log₂ κ(S) to log₂ κ(S) + 3 bits in every type
S9 marginalization	0.0	0.0	1.0	0.0	none: a gather
S10 online calibration	10.4	11.3	12.0, 5 flip	10.6	the window’s updates, as in S6e
C filter loop	9.5	10.6	12.0, 3/3 fail, 3 flip	9.4	the updates and the propagation, compounded over 54 operations (the run takes 46 min, nearly all of it posit16)

What the table shows:

• The losses are the problems’, not the types’. In every stage, float and posit32 lose about as many bits as double itself does (the long double column). The sensitivity is the conditioning of the problem the stage solves, not an algorithm that amplifies rounding. The total and stage bits agree to within a bit wherever a stage loses more than two, with two exceptions, both input rounding: S0 in posit32 (the pixel coordinates), and the noise-free S6a ladder, where the residual of the rounded input is not zero in double.
• The update is the precision-sensitive stage. S6e loses up to log₂ κ(S) + 3 bits. At σ = 0.05 px, κ(S) ≈ 2·10⁴, and that is 14 of a float’s 23 bits, and 14 of double’s 52. S10 inherits the same loss through its updates. A stage type narrower than float would land here first.
• Data movement is free. S3 (augmentation), S9 (marginalization) and S6c (QR compression) add no error of their own in any type, at any conditioning on the ladder: their stage bits lost are 0. posit16’s total of up to 1.8 bits is its rounding of the input.
• posit16 is too narrow for the estimator. It breaks S5, S6d and S6e at ordinary conditioning, and flips S10’s update decisions. It is usable only for the data-movement stages, and marginally for S0 and S1.
• posit32’s taper shows. It beats float wherever values are near 1, but loses more where they are not: S6b at 64 m, and the NIS in S6d.

These measurements choose the per-stage types for the end-to-end EuRoC check, the next step of #450.

Per-stage types on the estimator (#450)

The MSCKF backend stores its state in one type, T, and runs each stage in the type a stage-type policy names for it. The policy is a struct of type aliases, one per stage (sdk/msckf/stage_types.hpp). UniformStageTypes<T> puts every stage in T; it is the default, and compiles to the single-type backend. A mixed policy derives from it and overrides the stages it changes by name:

#include <branes/sdk/vio_estimator.hpp>

namespace ms = branes::sdk::msckf;

// Store the state in float, run the update (S5 → S6e) in double.
struct UpdateInDouble : ms::UniformStageTypes<float> {
    using update = double;
};

using Backend = branes::sdk::MsckfMixedBackend<float, UpdateInDouble>;
branes::sdk::VioEstimator<float, Backend> estimator(Backend(cameras));

Alias	Stage	What runs in that type
sensor_model	S0	pixel → normalized image point
propagation	S2	the mean and the covariance, per IMU sample
augmentation	S3	the clone and P’ = G P Gᵀ
update	S5 → S6e	triangulation, Jacobians, null space, compression, gate, EKF update
ekf_update (optional)	S6e	the EKF update alone; defaults to update
marginalization	S9	removal of the oldest clone
calibration	S10	the calibration states added to the fresh state

S1 (initialization) runs once, in the storage type, and S4 (the front end) hands the backend storage-type observations, so neither is in the policy.

ekf_update splits S6e — the Kalman gain, the mean correction and the covariance update — off the rest of the update. A policy that does not name it runs S6e in update, so using update = float; still puts the whole update in float. With it, S5 → S6d run in update and S6e in ekf_update, on the compressed measurement rounded into that type (s6_msckf_update::prepare is S5 → S6d without the update). If the innovation covariance that passed the gate in update does not factor in ekf_update, the track is rejected as an ill-conditioned innovation. S6e is the step a float update fails at on real data (below), which is why it has its own alias:

// The update in float, except S6e (gain, mean and covariance).
struct UpdateFloatEkfDouble : ms::UniformStageTypes<double> {
    using update = float;
    using ekf_update = double;
};

A stage in a type other than T runs on the state rebound into that type. The result is then rounded back to T (sdk/msckf/rebind.hpp). Each value is rounded to the nearest value of the target type. Rotations are renormalized in it, because a unit quaternion in one type is unit only to that type’s precision. A square-root covariance rebinds its factor, not P. The operators a stage consumes (the propagator, the camera updater, the camera model) are rebound once, when the backend is configured. Taps and the update observer still see the state in T.

Rebinding costs two copies of the state per stage call, O(n²) for P. S2 runs once per IMU sample, so a mixed propagation type costs about that much per sample. This is a diagnostic configuration, not a fast one.

Measured on the synthetic drive. tests/sdk/stage_types.cpp runs the estimator on an 8-frame synthetic drive with online calibration on. In each run one stage is in float and the storage and the rest are in double; the table gives the worst position deviation from the all-double run. For comparison, the whole estimator in float deviates by 1.4e-5 m on the same scene without calibration (tests/sdk/vio_estimator.cpp).

Stage in float	Worst |p − p_double|
S5 → S6e update	2.4e-5 m
S2 propagation	9.6e-6 m
S3 augmentation	2.8e-7 m
S0 sensor model	9.2e-8 m
S9 marginalization	1.9e-9 m
S10 calibration	3.1e-10 m

The update is the most sensitive stage here, as in the bench characterization above. The other stages do not rank as they do in the bench, because the bench measures one call and this run measures their effect on the trajectory:

• Propagation comes next because S2 runs once per IMU sample, so its rounding accumulates.
• S10 loses about 10 bits per call in the bench, but it runs once, when the state is seeded, so it barely moves the trajectory.
• S3 adds no error of its own in the bench, and S0 about 1.4 bits. What they cost here is mostly the rounding of their inputs into float: the whole state, for S3, on every clone; each observation, for S0.

Two more runs complete the picture. With every stage in long double and the state stored in double, the run deviates by 1.6e-14 m, double’s own rounding. With the state stored in float and the update in double, it deviates by 9.1e-6 m.

A state is only as PSD as its narrowest stage

Clone augmentation makes P exactly singular. A stage in float leaves those null directions at about −ε_float·‖P‖, even when the state is stored in double. Checking such a state against double’s PSD slack (16·ε·max diag) fails, although on the short synthetic drive the filter is fine. Check it at the ε of the narrowest type it passed through, as the test does. On a real run a float covariance update can do worse than this: see the EuRoC run below.

In the pipeline: --stage-types

vio_pipeline selects a policy by name with --stage-types NAME. The state is stored in double in every preset, so the front end, the output streams and the metrics are the same code in every run; only the stage arithmetic changes. The presets use native types only, because an emulated type end to end is about 10⁴× slower.

Preset	Stages in a type other than double
double	none (the default; identical to a run without the option)
float	every stage
long-double	every stage, in long double
sensor-model-float	S0
propagation-float	S2
augmentation-float	S3
update-float	S5 → S6e
marginalization-float	S9
calibration-float	S10, which runs only with online calibration (off in vio_pipeline)
float-update-double	every stage except S5 → S6e
ekf-update-float	S6e only
update-float-ekf-double	S5 → S6d
float-ekf-double	every stage except S6e

vio_pipeline --source euroc --dataset .../V1_01_easy/mav0 --stage-types update-float --out DIR

--stage-types applies to the MSCKF backend; with --invariant it is an error. So is a preset other than double with --capture or --live-assert: the stage recorder replays each boundary in double, so for a stage in another type it would record, or check, a run that did not happen.

First run, EuRoC V1_01. Every preset on the 2907 tracked frames of V1_01:

Preset	ATE	NIS / dof
double	0.2888 m	1.508
long-double	0.2888 m	1.508
sensor-model-float	0.2887 m	1.508
propagation-float	0.2844 m	1.508
augmentation-float	0.2671 m	1.507
marginalization-float	0.2882 m	1.508
float-update-double	0.2866 m	1.509
update-float-ekf-double	0.2626 m	1.528
float-ekf-double	0.2984 m	1.519
ekf-update-float	0.3452 m	1.500
update-float	128911 m (diverged)	0.476
float	128143 m (diverged)	0.037

Every stage but the update runs in float without harm: with all of them in float and the update in double, the trajectory moves by less than a centimetre of ATE. Don’t read the shifts of up to 2 cm in the other rows as a gain in accuracy: one run on one sequence can’t separate them from the run’s sensitivity to any perturbation. The update in float diverges.

Why the float update diverges on V1_01. A probe tapped every S6 call of the run, replayed it in double from the same input state, and swapped one sub-stage at a time into float. Three steps lead to the divergence:

1. Triangulation drops the near-static tracks. V1_01 starts with the vehicle almost still, and its first ~100 frames carry tracks with 0.03° of parallax. S5’s linear system A = Σ(I − d̂d̂ᵀ) then has a smallest eigenvalue of about the parallax angle squared, below float’s rounding, and its Cholesky factorization fails (135 of the 226 float failures; the rest fail in the Gauss-Newton refinement that starts from that solve). Over frames 0–110 the float run applies 577 updates where double applies 1196, and P grows to 1.67 where double’s is 0.19.
2. Takeoff collapses P in a few updates. At frame 105 one update cuts P’s largest diagonal from 1.49 to 0.135. The Joseph form (I − KH)P(I − KH)ᵀ + KRKᵀ forms the small posterior entries by cancellation of prior-sized terms, so its error is about ε·‖P_prior‖. In float that is larger than the posterior’s small eigenvalues: P comes out indefinite, its most negative eigenvalue −4·10⁻⁵ of its largest. Double, on the same updates, stays at its usual 10⁻¹¹.
3. An indefinite P runs away. At frame 111 a velocity variance goes negative (−1.8·10⁻⁶), and the next update multiplies it by 700. By frame 115 the negative eigenvalues are larger than the positive ones, and the estimate leaves the room.

Each step was checked by counterfactual:

Update run	ATE
all in float (update-float)	128911 m
S5 in double, the rest of the update in float	0.348 m
S6e (the EKF update) in double, the rest in float	0.256 m
all in float, parallax gate 0.1° (double with the gate: 0.394 m)	21674 m

From a healthy state a single float update is accurate. Replayed from the double run’s states, the float S6e update moves P by 6·10⁻⁷ relative (median; at most 10⁻⁴), and the float gate decides as double does on every track. The innovation covariance is well conditioned (κ ≈ 1.01: the measurement noise dominates it), so the gain is not the problem. What breaks is the covariance update on a large information gain. That is why the parallax gate does not help: it removes the float-only rejections, but it starves both runs, P drifts to 540 before takeoff, and the first collapse breaks the float update at once.

Two remedies follow. Keeping S6e, and with it the covariance update, in double is measured to work, and it is enough on its own: float-ekf-double runs every stage in float except S6e and stays at 0.298 m, and update-float-ekf-double (the counterfactual above, now a preset) at 0.263 m. S6e alone in float (ekf-update-float) does not diverge either, because double triangulation keeps P from growing before takeoff; it takes both halves to diverge. A square-root covariance, which is PSD by construction, keeps the whole update in float: with it, update-float runs V1_01 at 0.276 m (see the square-root form in float below).

Three sequences, with NEES

vio_trajectory_report takes the same --stage-types presets and adds the consistency metrics: RPE over 1 s, the 15-dof NEES (position and yaw anchored at the first post-init frame), and NIS. The runs below are on EuRoC V1_01 (easy, static start), MH_05 (difficult, static start) and V2_03 (difficult, moving start); NEES and NIS are per degree of freedom, and “diverged” means P went indefinite (NEES then skips those frames).

vio_trajectory_report --dataset .../MH_05_difficult/mav0 --label MH_05 --stage-types float-ekf-double --out DIR

Preset	V1_01 ATE · NEES · NIS	MH_05 ATE · NEES · NIS	V2_03 ATE · NEES · NIS
double	0.289 m · 436 · 1.51	0.755 m · 42.8 · 2.16	0.269 m · 140 · 14.7
long-double	0.289 m · 436 · 1.51	0.755 m · 42.8 · 2.16	0.269 m · 134 · 14.7
sensor-model-float	—	0.755 m · 42.7 · 2.16	0.275 m · 111 · 14.7
propagation-float	—	diverged	0.269 m · 111 · 14.7
augmentation-float	—	diverged	0.277 m · 110 · 14.9
marginalization-float	—	diverged	0.267 m · 187 · 14.9
update-float-ekf-double	0.263 m · 436 · 1.53	1.156 m · 72.5 · 2.02	0.303 m · 113 · 14.9
ekf-update-float	0.345 m · 446 · 1.51	diverged	0.281 m · 109 · 14.7
update-float	diverged	diverged	0.274 m · 111 · 14.9
float-ekf-double	0.298 m · 437 · 1.52	diverged	0.282 m · 109 · 14.9
float-update-double	0.287 m · 437 · 1.52	diverged	0.282 m · 131 · 14.9
float	diverged	diverged	0.369 m · 114 · 15.0

(V1_01’s single-stage rows other than the update are in the table above.)

The rule the three sequences show: a stage that writes P must not run in float. S2, S3, S9 and S6e write the covariance; S0 and S5 → S6d only read it. A stage in float runs on the state rounded into float, P included, and writes it back, which leaves an absolute error of about ε_float·‖P‖. That is harmless until the true small variances fall to that level, which happens when a burst of updates collapses P:

• MH_05 collapses P early: between frames 23 and 80 its largest diagonal drops 0.53 → 0.018 and its smallest (a gyro-bias variance) 4·10⁻⁴ → 9·10⁻⁷. With only propagation in float, the first negative variance appears at frame 78. Every stage that writes P fails here, each on its own, at the same point of the run.
• V1_01 collapses P less, at takeoff (frame 105), so each P-writing stage survives on its own; S6e fails only together with float triangulation, which lets P grow before takeoff.
• V2_03 starts moving, so there is no starved static phase and no collapse: every preset runs, the whole filter in float included (0.369 m).

The stages that only read P run in float on every sequence. With S5 → S6d in float (update-float-ekf-double) MH_05 costs 0.40 m of ATE (1.156 against 0.755 m) and its NEES rises from 42.8 to 72.5, so even the read-only part of the update is not free on the hardest sequence.

The NEES column also restates the known over-confidence (#212): 436, 43 and 140 per dof in double, where a consistent filter is near 1. Precision does not cause it: long double gives the same numbers on V1_01 and MH_05, and 134 against 140 on V2_03. How to keep a float fault from reaching the vehicle — the guards deployed autopilots run on P — is in Covariance Safeguards.

The square-root form in float (#533)

--covariance sqrt runs the square-root backend (#187), which carries a factor S with P = S·Sᵀ and updates it by QR instead of the Joseph form. It composes with every --stage-types preset.

vio_pipeline --source euroc --dataset .../MH_05_difficult/mav0 --covariance sqrt --stage-types float --out DIR

Until #533 the backend gave NaN on V1_01 even in double. The array update assumed a square S, but clone augmentation leaves S tall (81×75 on V1_01) until the next prediction, and that frame’s updates run in between. The update read past the end of each row of S and H·S, so S grew from frame to frame and, from about frame 108, went non-finite. Where it failed changed from run to run. With the pre-array sized by S’s columns, the square-root backend in double matches the full-covariance one on V1_01 to every printed digit.

Preset	V1_01 ATE · NIS	MH_05 ATE · NIS	V2_03 ATE · NIS
full double (reference)	0.289 m · 1.51	0.755 m · 2.16	0.269 m · 14.7
sqrt double	0.289 m · 1.51	0.755 m · 2.15	0.265 m · 14.9
sqrt ekf-update-float	0.288 m · 1.51	0.755 m · 2.15	0.282 m · 14.7
sqrt update-float	0.276 m · 1.54	0.761 m · 2.16	0.251 m · 15.0
sqrt float-ekf-double	0.425 m · 1.50	0.768 m · 2.17	0.303 m · 14.9
sqrt float	0.295 m · 1.51	0.784 m · 2.22	0.282 m · 15.0

Each of the five presets measured holds on these three sequences, the whole filter in float included. The full-covariance form diverges with update-float on V1_01 and MH_05, and with any P-writing stage in float on MH_05. So on these runs the rule above was about P, not about the stages: a stage that writes the covariance needs double when the covariance is P, and float was enough when it was S. Three runs show that much and no more; they are not a guarantee for other sequences or other float stages.

Why the factor does better: it carries the covariance in half the dynamic range (κ(S) = √κ(P)), and the update never forms P, so there is no cancellation of prior-sized terms like the Joseph form’s. S·Sᵀ is positive semidefinite in exact arithmetic. Computed in floating point it can still have small negative eigenvalues from rounding, and S itself can become non-finite (as #533 did), so the health check the safeguards report recommends applies to the square-root form too. In these runs the takeoff collapse that breaks the float Joseph form left the factor intact.

MH_05 is the hard case: all in float costs 3 cm of ATE (0.784 against 0.755 m). The V1_01 row for float-ekf-double, 0.425 m, is the one outlier: the S6e array update in double runs on a factor that the float stages wrote. One run per preset can’t tell whether that is a precision effect or the run’s sensitivity to any perturbation (V2_03 in double moves 4 mm between the two forms, which differ only in rounding).

The MTL5 gap register (#451)

Every operator the pipeline needs that MTL5 doesn’t provide is recorded with its pipeline sites, the interim workaround, and an upstream issue (stillwater-sc/mtl5). Known so far:

• scalar × matrix for non-arithmetic scalars (mtl5#538);
• a vectorizable generic dense product (mtl5#549) and a CRS × dense-matrix product (mtl5#550), both in mtl5#552: done, MTL5 5.14;
• a non-owning CRS view over caller arrays (compressed2D copies them);
• block insert/scatter into dense2D;
• a symmetric A·P·Aᵀ congruence (a gemmt-style triangular-output product; syrk/symm exist but don’t fit, see above) (mtl5#656);
• a left-null-space helper;
• an opt-in term-by-term ordering in cholesky_factor, and NaN-pivot rejection (mtl5#556): done, MTL5 5.14;
• an opt-in textbook Householder reflector for QR and the null space (mtl5#559): done, MTL5 5.15. Multi-RHS solves need nothing new: two trsm calls do it;
• a matrix exponential.

Lie-group exp/log and Jacobians stay in branes::math::lie; they are not MTL5’s job.