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 theblas_traitstriple-product block callverifyProductwith the magnitudes of the intermediates —|m_1|\,|m_1|^Tand|\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, soS\,(SS)^Tis bounded by|S|\,(|S|^2)^T— for a nilpotentSthe two differ in kind,|S|^3 = 0while the two evaluation orders still disagree by a rounding error.test/product_test_helpers.h:product_toleranceloses its last two callers and is removed.product_error_boundnow returnsdoubleandverifyProductaccumulates the error norm indouble, with\varepsilonstill the operand scalar's: narrowed tohalf, the triple-product bound is already infinite atn = 160with the test's own random operands, andinf <= infthen accepted any result. A non-finite bound or error now fails.test/product_small.cpp: deterministic regressions for both — the nilpotent3\times 3case, and an all-inf/all-NaNhalfresult against a bound that overflowshalf(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.0685Failure 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