[RobMina]

Summary

Every conversion between the parameter (covariance) and weight (inverse covariance) form of a fit state goes through
FitData::invert(), which inverts the symmetric 6x6 matrix with SMatrix::Invert(), i.e. Bunch-Kaufman diagonal
pivoting. In Mu2e reconstruction this is now the single most expensive function wherever a KinKal fit runs
(callgrind, Offline main + KinKal main with #227):

Workload (Mu2e reco)	Inversion share of all event-processing instructions	Calls / event
on-spill CE + 1BB pileup, LoopHelix fits	14.9 %	43,600
off-spill cosmics, LoopHelix + CentralHelix fits	17.7 %	65,400
extracted-detector cosmics, KinematicLine	14.1 %	2,400

Two thirds of the weight-side conversions come from FitState::wData() in DomainWall::process and
Material::process.

The matrices are covariances or their inverses, so in a healthy fit they are symmetric positive-definite. This PR
inverts them with SMatrix::InvertChol() (ROOT CholInverter<6> / CholeskyDecomp<double,6>, the fixed-size
decomposition) and falls back to the existing Invert() on the original matrix when the decomposition reports a
non-positive pivot:

void invert() {
  // Cholesky for the symmetric positive-definite case; Bunch-Kaufman on the original matrix if it is not
  DMAT orig(mat_);
  if(!mat_.InvertChol()){
    mat_ = orig;
    if(!mat_.Invert()) throw std::runtime_error("Inversion failure");
  }
  vec_ = mat_*vec_;
  // finiteness checks unchanged

So a matrix that Bunch-Kaufman inverts today is still inverted, and one it cannot invert still throws
"Inversion failure". The other Invert() calls in KinKal (Parameters::delta, ParameterHit, Track::initFitState)
are not changed; they are rare.

Why the fallback is needed

On 30,000 matrices dumped from FitData::invert in Mu2e on-spill reconstruction (first 3 events, 15,165 distinct),
Cholesky refuses 42 (0.14 %). All 42 have exactly one negative eigenvalue (from -0.1 to -1.3e5, the others 0.07 to
1.4e6), and in 9 of them a diagonal element is itself negative, so they are not proper covariance or weight
matrices; Bunch-Kaufman inverts them without complaint today. A likely source is Hit::unbiasedParameters(), which
subtracts the hit's weight from the reference weight; I have not confirmed which caller produces them. In full reco
the fallback is taken for 0.08 % (on-spill) and 0.10 % (off-spill) of calls.

Tests

• Standalone, header-only (ROOT 6.38 SMatrix, g++ 13.3 -O3):

  Input	Cholesky refused	Relative difference of the two inverses, p50 / p99 / max	Cost per call (copy + invert)
  30,000 real matrices from Mu2e reco	42 (0.14 %)	3.6e-14 / 1.4e-8 / 3.0e-7	207-210 -> 110-111 ns (0.53x)
  synthetic SPD, condition 1e4 / 1e8 / 1e12	0	max 1.1e-12 / 9.2e-9 / 1.1e-4	239-242 -> 105-106 ns (0.44x)
  callgrind, 10,000 SPD matrices			3,361 -> 753 instructions (0.22x)

  On badly conditioned synthetic matrices Cholesky's backward error |A·A⁻¹ − I| is larger than Bunch-Kaufman's (2x at
  the median, up to 60x at the tail for condition 1e12); on the real matrices the two agree to ~1e-14 typically.

• KinKal unit tests: 26/26 pass.

• Mu2e reconstruction, cost (callgrind, same 400 events before and after; base = unpatched build of the same
  commits):

  	on-spill CE + 1BB	off-spill cosmics
  instructions per event	928 M -> 821 M (-11.5 %)	1,169 M -> 1,003 M (-14.2 %)
  LoopHelix fit modules	-24 %	-22 %
  FitData::invert, instructions per call	3,308 -> 811	~3,300 -> 810
  inversion share of the job	15.5 % -> 4.3 %	17.7 % -> 5.2 %

  Same host, same 2000 on-spill CE + 1BB events: each LoopHelix fit module is 6 % faster on the mean and 7 % on the
  median; the median event is 3.4 % faster, the mean event 16 % (the host was shared, so the mean is the noisier
  number); output, calorimeter and selection modules are unchanged. The time gain is smaller than the instruction
  gain, as in the standalone timing.

• Mu2e reconstruction, physics output. The arithmetic changes, so fits are not bitwise identical; the comparison is
  statistical, same events in both arms:

  Sample	Events	Tracks before / after	Truth-matched efficiency ratio (after/before)	Largest shift of a per-track mean
  on-spill CE + 1BB, 4 LoopHelix hypotheses	24,950	99,000 / 99,018	1.0002 (+1.0 sigma)	1.0 sigma
  off-spill cosmics, CentralHelix	10,000	1,241 / 1,241	1.0000 (0 sigma)	1.1 sigma
  extracted-detector cosmics, KinematicLine	28,795	20,758 / 20,757	0.9999 (-1.0 sigma)	0.2 sigma

  Per-track means compared: fit status, fit consistency, chi2, ndof, active hits, material crossings, trajectory
  segments. On tracks in events with the same track count, chi2 and fit consistency change in the last digits for
  89-100 % of tracks (CentralHelix fit consistency 89 %, all others 96-100 %), as expected. ndof is unchanged for
  KinematicLine. ndof changes for 0.6 % of CentralHelix tracks and 29 % of LoopHelix tracks: the
  Mu2e LoopHelix fit's drift/background hit-state updates (ANN) can flip a hit for a tiny change of the fit, and the
  flips are balanced (658 of 24,950 events change their track count: 317 with more tracks before, 341 after). No
  exceptions in either arm.

Not changed

• Parameters::delta, ParameterHit and Track::initFitState still call Invert() directly.
• No change to which matrices are produced. Whether the indefinite ones (0.1 % of conversions) should exist at all,
  and whether DomainWall::process could stay in weight space and save one of its two inversions per wall, are
  separate questions.

[oksuzian]

PR Review Summary — #229

Reviewed at head a9a8aff03f02037ecf08fabe8c18bd5f4f98547d. first pass

Decision

• 🟢 approve. I reproduced the speed-up, and on the KinKal toy fits the results are statistically identical to main. One change in failure behaviour is worth a comment or a guard (finding 1). It does not block, because I could not find a way for a normal KinKal fit to reach it.

Scope understood

• FitData::invert() now tries SMatrix::InvertChol() first. If Cholesky refuses, it inverts the original matrix with Invert() (Bunch-Kaufman), as before. The finiteness check after the inversion is unchanged.
• One file, General/FitData.hh. The other direct Invert() calls, in Parameters::delta, ParameterHit::chisq and Track::setStatus, are unchanged.

Findings

1. 🟡 [S2] An exactly singular, positive semi-definite matrix no longer throws.

   • Evidence: in ROOT 6.38 CholeskyDecomp<double,6> (Math/CholeskyDecomp.h, _decomposer<F,6,M>), the only test is pivot <= 0 at each step. Rounding often leaves the last pivot of a singular PSD matrix a tiny positive number, so the decomposition reports success.
   • Measured, standalone:
     • 20,000 random rank-5 PSD 6x6 matrices: Bunch-Kaufman refuses 1.6% and Cholesky refuses 55%.
     • But 124 of them are accepted by Cholesky and refused by Bunch-Kaufman. For those, the "inverse" has elements up to 8.5e16 and |A·A⁻¹ − I| up to 86.
     • Rank-4: 17 such cases.
     • Through FitData::invert() itself, a matrix with row 5 equal to row 4 throws Inversion failure on main, but on this PR returns a finite matrix with elements around 1e15 and |A·A⁻¹ − I| = 1.
   • Impact: as far as I can tell, nothing reachable today.
     • The fit state always starts from the de-weighted seed covariance (Track::initFitState, dwt_ times the seed), and everything after that adds weight. So the states that DomainWall/Material::process invert are full rank.
     • In the toy fits below, Cholesky refused only 43 of 7.6 million inversions, and nothing hit this case.
     • It would matter for a caller that hands FitData a singular covariance, for example a seed with a zero-variance parameter. That used to throw; now a huge finite inverse goes into the fit without any message. The wiki entry for issue #225 W2 recorded "an exactly singular matrix makes Invert() return false" as the guard here.
   • Suggested fix, any one of:
     • Keep the old behaviour by also rejecting a Cholesky result whose smallest pivot is below a relative tolerance. CholeskyDecomp exposes the factor, so this costs a few compares.
     • Check |diag(A)·diag(A⁻¹)| against a bound.
     • Or document in the comment that singular PSD input is no longer an error.

2. ⚪ [S3] Cholesky is less accurate than Bunch-Kaufman on badly conditioned SPD matrices, as your description says. Measured on synthetic matrices with KinKal-like units (phi0 scaled by 1e-3):

   condition number	relative difference of the two inverses, p50 / max	backward error, Cholesky p50 / max	backward error, Bunch-Kaufman p50 / max
   1e8	5.0e-10 / 9.5e-9	9.5e-7 / 2.4e-4	4.8e-7 / 9.4e-6
   1e12	5.7e-6 / 9.9e-5	1.3e-2 / 13	5.9e-3 / 0.078
   1e16	6.1e-2 / 24	1.3e2 / 2.7e5	64 / 1.0e4

   The median backward error is about 2x worse, and the tail is up to 25-170x worse. At realistic conditioning this does not show up in the fit (A/B below). Worth one line in the comment, so nobody later concludes the two are interchangeable for ill-conditioned input.

🟢 Verified by running code (ROOT 6.38.00, g++ 13.3, -O3, KinKal main 14fd223 vs this head, both built with cmake Release):

• NaN and Inf still throw. Cholesky reports success on a matrix with a NaN element: NaN <= 0 is false, so the pivot test passes. I measured 2000/2000 cases with non-finite output. But the unchanged finiteness check then throws Inversion failure in both arms. One Inf off the diagonal is refused by Cholesky in 3% of cases and otherwise produces NaN, and it throws in both arms. An all-zero matrix throws in both arms.

• Indefinite matrices take the fallback, every time. A single negative eigenvalue anywhere from -1e-12 to -1, or a negative diagonal element, is refused by Cholesky in 20,000 of 20,000 cases. Bunch-Kaufman then inverts it, and the result is bit-identical to main through FitData::invert().

  • This is the case your 42 refused Mu2e matrices fall in.
  • It is what Hit::unbiasedParameters() (Detector/Hit.hh:61-70, reference weight minus the hit's weight) produces when it over-subtracts. A simulated over-subtraction of 1.000001x to 2x the information along a hit direction gives 100% refusal and a clean fallback.
  • At exactly 1x (singular), Cholesky accepts 45% of cases; see finding 1.

• Well-conditioned SPD matrices agree to rounding. At condition about 1e2, the relative difference is 4.3e-16 at the median and 7.5e-15 at most, with the same backward error.

• Speed. Over 3 runs, copy plus InvertChol (this PR) takes 141-158 ns per call against 313-342 ns for Invert() on main, a ratio of 0.45-0.46. The extra DMAT orig(mat_) copy costs 8 ns, about 5% of the new cost, so it is negligible. This matches your 0.44-0.53x.

• KinKal unit tests: 26/26 pass at both main and this head.

• Fit A/B, from the KinKal toy drivers: same seed 4242, 2000 events each, the drivers' own default schedules (MCAmbigFixedField for the helices, driftfit with a momentum constraint for KinematicLine).

  	LoopHelix	CentralHelix	KinematicLine
  fit status identical	2000/2000	2000/2000	2000/2000
  ndof identical	2000/2000	2000/2000	1996/2000
  iterations identical	2000/2000	2000/2000	2000/2000
  signed Δχ² mean (converged)	+0.0001 ± 0.0006	−0.0001 ± 0.0002	−0.0034 ± 0.0053
  χ²/ndof mean, main vs PR	0.9858 / 0.9858	0.9869 / 0.9869	0.9430 / 0.9430
  momentum pull mean, main vs PR	−0.1466 / −0.1473	−0.1774 / −0.1775	−0.1080 / −0.1040
  largest |Δχ²|	0.75	0.22	8.5, when one hit changes state (ndof 72→70)

  The fits are not bit-identical, as you said, but there is no shift in χ², χ²/ndof or momentum pull. The per-event spread in KinematicLine comes from 4 fits where a drift-ambiguity hit flips. KinematicLine ends with status −3 (pull tolerance) in both arms, so that is not from this PR.

• The fallback is rare in the toy fits. With review-only counters in a copy of this head: LoopHelix 0 of 2.69M inversions refused, CentralHelix 0 of 2.38M, KinematicLine 43 of 7.63M (0.0006%). Your Mu2e numbers, 0.08-0.10%, are higher, which fits the extra effects in Mu2e reco.

Validation check

• Build/tests run: GitHub Actions Tests pass on all 4 configurations at this head. Locally: both arms built, ctest, the standalone tests and the fit A/B as above. Scripts are in claude-scratch/kinkal229_review/ (run_inv_test.sh, run_edge.sh, run_singular.sh, run_ab.sh, run_compare.sh, run_count.sh).
• Not run: Mu2e reconstruction with this change. The Mu2e cost and physics tables in the description are yours.
• Config contract check: n/a, no configuration change.
• Cross-repo consistency: no interface change. Offline picks this up with the next KinKal tag and envset.

Residual risk

• Finding 1, if a future caller passes a singular covariance.
• The accuracy loss in finding 2 at condition numbers of 1e12 and above, if KinKal fits ever operate there. I did not measure the condition numbers of real Mu2e fit states.

Author follow-ups

1. Optional: a pivot tolerance, or a comment, for the singular PSD case (finding 1).
2. Optional: one line in the comment on the accuracy trade-off (finding 2).

[RobMina]

I expanded the doc comment as requested. The problem case identified in the AI review did not come up in testing and seems extremely unlikely.

[oksuzian]

PR Review Summary — #229

Reviewed at head 5a8f0c73a32c8818dcb71397735132f2b6038d45. re-review of a9a8aff

Decision

• 🟢 approve

Scope understood

• One new commit, 5a8f0c73, expands the comment above FitData::invert() in General/FitData.hh. No code changed since my approval at a9a8aff0.

Prior findings

• 🟢 [was S2] An exactly singular positive semi-definite matrix no longer throws — addressed in 5a8f0c73 by the documentation option I listed. The comment now says such a matrix can pass Cholesky and return a huge finite inverse, and that fit states are full rank, so the case does not arise in a fit. That matches what I measured: no occurrence in 7.6 million toy-fit inversions. A future caller that hands FitData a singular covariance now has a stated warning rather than a silent change.
• 🟢 [was S3] Cholesky is less accurate than Bunch-Kaufman on badly conditioned input — addressed in the same comment, which now states the trade-off.

Validation check

• Build/tests run: none at this head, since only a comment changed. The repository's four test jobs (Debug and Release, Linux and macOS) are green at this head.
• Config contract check: not applicable.
• Cross-repo consistency: no Mu2e change needed.

Residual risk

• Unchanged from my first pass: none reachable in current Mu2e fits.

[brownd1978]

Thanks for this improvement Rob!