Tests: Bound the product associativity checks by |A|*|B| instead of the result

product_large and product_small flake at roughly one failure per 400 calls to product() because the associativity and blas_traits blocks in test/product.h measure their error relative to the result, while both compare quantities that cancel. Moving those two sites onto the absolute Higham–Mary bound that product_small and product_trmv already use — verifyProduct(actual, expected, |A|, |B|, num_products) — takes the measured failure count from 53 in 200 000 draws to 0 in 2.6 million checks across six scalar types.

On master, 3 of 600 product_large_2 r10 s<n> runs fail (seeds 42, 182, 184); all three pass on this branch. Most recently seen in CI as job 16357381683 (aarch64, gcc-10, seed 1937694149):

Difference too large wrt tolerance 2.9976e-14, relative error is: 1.37439e-13
Test product(MatrixXd(internal::random<int>(1, 10), internal::random<int>(1, 10))) failed in ../test/product.h (120)
    verifyIsApprox((m1 * m1.transpose()) * m2, m1 * (m1.transpose() * m2), tol)

Why no relative tolerance can work here

Both sides of the associativity check chain two products, so by Higham & Mary each carries an error of order

\lambda\sqrt{k}\,\varepsilon\,\bigl\||m_1|\,|m_1|^T|m_2|\bigr\|_F ,

while the result itself vanishes with m_1^Tm_2. The amplification

\frac{\bigl\||m_1|\,|m_1|^T|m_2|\bigr\|_F}{\bigl\|m_1m_1^Tm_2\bigr\|_F}

is therefore unbounded; 1.3\times10^{4} was observed in 200 000 draws of the shape the test uses. product_tolerance is a relative model, so no choice of its constant fixes these two sites — the bound has to be taken against the magnitude of the intermediates. verifyIsApproxScaled in test/numerical_test_helpers.h already records the same argument in prose.

Measured over 200 000 draws of product(MatrixXd(random(1,10), random(1,10))):

bound failures worst error/bound
product_tolerance, i.e. master 53 / 200 000 (0.027 %) 48
verifyProduct against $`| m_1

What changed

  • test/product.h: the associativity block and the blas_traits triple-product block call verifyProduct with the magnitudes of the intermediates — |m_1|\,|m_1|^T and |\texttt{square}|^2 — as the bound's operands. The scaled variant passes |s_1|\,|\texttt{square}| so the bound carries the scale factor. Conjugation leaves the magnitude chain at |S|^3; transposition preserves each entry's magnitude but reorders the chain, so S\,(SS)^T is bounded by |S|\,(|S|^2)^T — for a nilpotent S the two differ in kind, |S|^3 = 0 while the two evaluation orders still disagree by a rounding error.
  • test/product_test_helpers.h: product_tolerance loses its last two callers and is removed. product_error_bound now returns double and verifyProduct accumulates the error norm in double, with \varepsilon still the operand scalar's: narrowed to half, the triple-product bound is already infinite at n = 160 with the test's own random operands, and inf <= inf then accepted any result. A non-finite bound or error now fails.
  • test/product_small.cpp: deterministic regressions for both — the nilpotent 3\times 3 case, and an all-inf/all-NaN half result against a bound that overflows half (fails at the parent commit).

The constants are unchanged from what the call sites asked for: num_products 3 and 4, lambda 5, so the bounds are 15\sqrt{k}\,\varepsilon and 20\sqrt{k}\,\varepsilon times the magnitude of the intermediates. Against the worst observed ratio of 2.6 that is a margin of 5.8, comparable to the margin product_small runs with. With the corrected chain the transposed predicate's worst ratio is 1.39 over 470 000 draws (zero failures), and half at n \le 320 — previously vacuous — is now checked at a worst ratio of 0.12.

Validation

product_large_{1..5,7..12}, product_small_{1..6,8,10,50..55} and product_trmv_{1..4} pass at r20. Seeds 42, 182 and 184 — the 3 failures in a 600-seed sweep of product_large_2 r10 at the parent commit c93cc5c0b — all pass here. The sweeps below are standalone programs built against this branch's headers, not CI jobs.

Appendix: calibration and stress sweeps

Host: 13th Gen Intel(R) Core(TM) i7-13700HX under WSL2, g++ (Ubuntu 13.3.0-6ubuntu2~24.04.1) 13.3.0 -O2.

Worst observed error / (eps * sqrt(k) * || |A| |B| ||_F) per scalar type and size range; the shipped bounds multiply that denominator by 15 (associativity) and 20 (blas_traits):

double           maxdim= 10  site1 max=1.88   site2 max=2.6
float            maxdim= 10  site1 max=2.01   site2 max=2.49
complex<double>  maxdim= 10  site1 max=1.6    site2 max=1.95
complex<float>   maxdim= 10  site1 max=1.68   site2 max=1.94
bfloat16         maxdim= 10  site1 max=1.79   site2 max=1.38
half             maxdim= 10  site1 max=1.52   site2 max=1.27
double           maxdim= 20  site1 max=1.74   site2 max=1.19
float            maxdim= 20  site1 max=1.71   site2 max=1.47
complex<double>  maxdim= 20  site1 max=1.58   site2 max=1.74
complex<float>   maxdim= 20  site1 max=1.23   site2 max=2.05
bfloat16         maxdim= 20  site1 max=1.87   site2 max=1.06
half             maxdim= 20  site1 max=1.51   site2 max=1.27
double           maxdim=320  site1 max=0.427  site2 max=0.224
float            maxdim=320  site1 max=0.441  site2 max=0.521
complex<double>  maxdim=320  site1 max=0.861  site2 max=0.911
complex<float>   maxdim=320  site1 max=0.179  site2 max=0.571
bfloat16         maxdim=320  site1 max=0.301  site2 max=0.0259
half             maxdim=320  site1 max=0      site2 max=0.0685

Failure counts for the five verifyProduct predicates exactly as this branch spells them, 520 000 draws (2.6 M checks):

double           maxdim= 10  N=100000  fails=0
float            maxdim= 10  N=100000  fails=0
complex<double>  maxdim= 10  N=25000   fails=0
complex<float>   maxdim= 10  N=25000   fails=0
bfloat16         maxdim= 10  N=5000    fails=0
half             maxdim= 10  N=5000    fails=0
double           maxdim= 20  N=100000  fails=0
float            maxdim= 20  N=100000  fails=0
complex<double>  maxdim= 20  N=25000   fails=0
complex<float>   maxdim= 20  N=25000   fails=0
bfloat16         maxdim= 20  N=5000    fails=0
half             maxdim= 20  N=5000    fails=0
TOTAL FAILURES: 0
Edited by Rasmus Munk Larsen

Merge request reports

Loading
Loading