Core: Refactor approximate comparisons and fix false positives
This makes dense isApprox and isMuchSmallerThan correct at extreme scales and under flush-to-zero (FTZ), and fixes the shared test assertions that accepted opposite infinities or dropped NaNs. On master, a.isApprox(2*a) is true for nonzero a at scales 1e200 and 1e-200 because the squared norms overflow or underflow. Comparisons of 256 or more coefficients run at the same speed; small fixed-size ones cost up to about 1.5 ns more.
Core
- Ordinary path. Squared norms are compared as before, but only when nothing overflows and
\min(\|x\|^2, \|y\|^2) \ge n\,\lambda / \varepsilon_{\text{acc}}, where\lambdais the smallest normal number and\varepsilon_{\text{acc}}is the accumulator's epsilon. Under that condition, squares that underflow or are flushed change the result by at most\varepsilon_{\text{acc}}relative. Half and bfloat16 are widened to float before squaring, so half comparisons at the default precision stay on this path. Callgrind counts 4.4–10× fewer instructions per comparison for half and 2.5–8.4× fewer for bfloat16 (16–4096 coefficients, GCC 13-O3, Intel Core i7-13700HX). Float and double code is unchanged. - Scaled fallback. All other inputs go to an out-of-line
SafeScalingpath, so the ordinary path still inlines. Operands whose largest component satisfiesM < 1are scaled up by a power of two before the subtraction, so under FTZ it loses only differences belowM\,\lambda. On master,Vector2f(1.5f * FLT_MIN, 0).isApprox(Vector2f(FLT_MIN, 0), 0)is true under FTZ, on the host and with CUDA-ftz=true. NaNs propagate. Operands that are not dense-shaped, such asHomogeneous, are evaluated first. Dense expressions stay lazy, so comparisons do not allocate. - Custom scalars keep the existing algebraic path. Comparing half with float still fails to compile; the new
approx_mixed_scalarsfailtest checks this. Fuzzy.handSparseFuzzy.hare renamed to the private headersApproxComparisons.handSparseApproxComparisons.h. Sparse behavior is unchanged.dummy_precision()is now documented as a default relative tolerance, not an error bound.functor_traits<scalar_unary_pow_op>no longer advertises packet access for a complex exponent. Only the real-exponent specializations definepacketOp, soArrayXcf r = a.pow(std::complex<float>(2, 1))withArrayXf afails to compile on master in vectorized builds.
Public names and default tolerances are unchanged. Two results differ from master on purpose: x.isMuchSmallerThan(x) and x.isApprox(y, +inf) return false when x contains an infinity. Master returns true for both because it evaluates \infty \le \infty.
Tests
- The scalar, packet, matrix and array comparison helpers now live in
numerical_test_helpers.h. They reject opposite infinities, propagate NaNs, and give correct relative-error diagnostics for integers, zero inputs and custom complex scalars. Packet tests keep theirdummy_precision()default. test_isApproxWithRefis removed: it compared shifted values, which could round the difference away.adjointnow uses exact integer checks or explicit bounds built from epsilon, dimension and operand magnitude. The 2×2 extreme-scale Bunch–Kaufman tolerance tightens from1e-3to8 * epsilon. The structured-matrix tests share one IEEE reference product and one comparator.- Some checks passed only because both squared norms underflowed to zero:
adjoint's zero-norm checks now evaluate the subnormal vector once. ARMv7 NEON flushes subnormal products evaluated in packets but keeps them in scalar reads of the same lazy expression.- The SVD
check_solverbaseis skipped whenrank()'s clamp atnumeric_limits::min()drops a positive singular value\gethreshold() * sigma_0, becausem \cdot \mathrm{solve}(m x) = m xcannot hold then.svd_least_squarecovers those inputs. Zero matrices keep the check (jacobisvd_56,bdcsvd_45).
cwise_min_maxworks around a GCC 10 AArch64 miscompile: the compiler kept the unreduced 128-bit partial result ofmaxCoeff()inv8across a call, but AAPCS64 preserves only the low 64 bits.- New:
approx_comparisons(four parts, covering FTZ, extreme scales,Homogeneous, lazy sums, complex-exponentpow,CustomComplex, and the no-allocation property),numerical_test_helpers, andbenchmarks/Core/bench_approx_comparisons.cpp.
Thanks to @onalante-ebay and @florian360 for their reviews.
Performance
Case (Apple M4, Clang 17 -O3, NEON) |
Master ns | This MR ns | Ratio |
|---|---|---|---|
isApprox<double>, fixed 4 coefficients |
1.063 ± 0.002 | 1.378 ± 0.002 | 1.30× |
isMuchSmallerThan<double>, matrix reference, 16 coefficients |
2.758 ± 0.009 | 4.246 ± 0.020 | 1.54× |
isApprox<double>, 4096 coefficients |
2008.208 ± 5.512 | 2007.347 ± 4.451 | 1.00× |
Values are the median ± MAD of ten CPU-time samples (two alternating runs of five repetitions each), with affinity and clocks not locked. Operands are x = ones and y = (17/16) * ones with prec = 1/8; the smaller-than cases use x = ones/32. These timings come from the branch's first revision. Later revisions change only the scaled fallback, the half and bfloat16 paths, and the classification of non-finite values in the guard.
All Apple M4 measurements
| Comparison | Coefficients | Master ns | This MR ns | Ratio |
|---|---|---|---|---|
isApprox, float |
4 (fixed) | 0.814 ± 0.001 | 1.127 ± 0.003 | 1.38× |
isApprox, float |
16 | 3.992 ± 0.022 | 4.123 ± 0.104 | 1.03× |
isApprox, float |
256 | 37.761 ± 0.178 | 37.693 ± 0.172 | 1.00× |
isApprox, float |
4096 | 970.167 ± 1.767 | 971.549 ± 1.862 | 1.00× |
isApprox, double |
4 (fixed) | 1.063 ± 0.002 | 1.378 ± 0.002 | 1.30× |
isApprox, double |
16 | 5.262 ± 0.008 | 5.773 ± 0.010 | 1.10× |
isApprox, double |
256 | 80.805 ± 0.230 | 80.781 ± 0.162 | 1.00× |
isApprox, double |
4096 | 2008.208 ± 5.512 | 2007.347 ± 4.451 | 1.00× |
isApprox, std::complex<double> |
4 (fixed) | 2.757 ± 0.008 | 3.139 ± 0.003 | 1.14× |
isApprox, std::complex<double> |
16 | 12.072 ± 0.030 | 12.410 ± 0.022 | 1.03× |
isApprox, std::complex<double> |
256 | 176.557 ± 0.292 | 177.047 ± 0.311 | 1.00× |
isApprox, std::complex<double> |
4096 | 2821.382 ± 5.560 | 2823.367 ± 4.781 | 1.00× |
isMuchSmallerThan (matrix reference), double |
4 (fixed) | 0.626 ± 0.001 | 0.940 ± 0.002 | 1.50× |
isMuchSmallerThan (matrix reference), double |
16 | 2.758 ± 0.009 | 4.246 ± 0.020 | 1.54× |
isMuchSmallerThan (matrix reference), double |
256 | 45.239 ± 0.126 | 46.603 ± 0.315 | 1.03× |
isMuchSmallerThan (matrix reference), double |
4096 | 1289.924 ± 4.226 | 1287.862 ± 1.456 | 1.00× |
isMuchSmallerThan (scalar reference), double |
4 (fixed) | 0.398 ± 0.000 | 0.571 ± 0.004 | 1.44× |
isMuchSmallerThan (scalar reference), double |
16 | 1.497 ± 0.012 | 2.015 ± 0.011 | 1.35× |
isMuchSmallerThan (scalar reference), double |
256 | 22.743 ± 0.075 | 23.124 ± 0.050 | 1.02× |
isMuchSmallerThan (scalar reference), double |
4096 | 643.944 ± 1.558 | 645.755 ± 1.532 | 1.00× |
Local validation
| Configuration | Result |
|---|---|
GCC 13.3.0, C++14, -O3, x86-64 |
approx_comparisons_1–4, numerical_test_helpers, adjoint_1–24, mixingtypes, geo_homogeneous_1–3, bunchkaufman_1–7, eigensolver_complex_6, array_cwise_19–21 passed |
| Apple clang 21, C++14, Release, NEON | approx_comparisons_1–4, numerical_test_helpers, jacobisvd_11/13/56, bdcsvd_7/16/45, complex_sqrt, eigensolver_complex_6, vectorwise_norm_12, matrix_exponential_3 passed; seed 19, 20 repetitions |
| LLVM Clang 17, C++14, Release / Debug with vectorization disabled and row-major default | 66 / 7 CTest cases passed, covering adjoint, arrays, geometry, sparse, structured, Bunch–Kaufman, packet math (real, complex, half, bfloat16) and matrix functions; seed 19, 20 repetitions |
| Apple Clang, ASan + UBSan | Comparison and helper regressions passed |
CUDA 12.8 sm_89 and CUDA 13.0 |
Device isApprox for float, double, complex and half at scales 1, 2^{\text{min\_exponent}+4} and 2^{\text{max\_exponent}-8} passed with and without -ftz=true; gpu_basic passed |
-O3 -ffast-math |
Finite extreme-scale checks passed. This MR makes no IEEE special-value claim under that flag |
| Master | The extreme-scale dense regression, the float FTZ case, the complex-exponent pow assignment, and the opposite-infinity and NaN helper probes fail |
| Formatting and static checks | clang-format 17, git diff --check, repository style checks and clang-tidy passed |
The GCC 13.3 and CUDA 12.8 rows ran on the final revision; the others ran on earlier revisions of this branch. HIP, SYCL and the Doxygen build were not run locally.
Known issues not fixed here
MatrixBase::isOrthogonalreports a nonzero vector as orthogonal to itself at scales1e100and1e-100.- Sparse
isApproxand complex scalarinternal::isApproxaccept twofold differences at scales1e200and1e-200. - Real scalar
internal::isApprox, whichisApproxToConstantalso uses, accepts opposite infinities. This MR guards only the test wrappers. - Other tests use loose residual tolerances, take maxima that can hide NaNs, or downgrade numerical assertions to warnings. Tightening each one needs its own error model and an independent reference.