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

  1. 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 \lambda is 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.
  2. Scaled fallback. All other inputs go to an out-of-line SafeScaling path, so the ordinary path still inlines. Operands whose largest component satisfies M < 1 are scaled up by a power of two before the subtraction, so under FTZ it loses only differences below M\,\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 as Homogeneous, are evaluated first. Dense expressions stay lazy, so comparisons do not allocate.
  3. Custom scalars keep the existing algebraic path. Comparing half with float still fails to compile; the new approx_mixed_scalars failtest checks this.
  4. Fuzzy.h and SparseFuzzy.h are renamed to the private headers ApproxComparisons.h and SparseApproxComparisons.h. Sparse behavior is unchanged. dummy_precision() is now documented as a default relative tolerance, not an error bound.
  5. functor_traits<scalar_unary_pow_op> no longer advertises packet access for a complex exponent. Only the real-exponent specializations define packetOp, so ArrayXcf r = a.pow(std::complex<float>(2, 1)) with ArrayXf a fails 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 their dummy_precision() default.
  • test_isApproxWithRef is removed: it compared shifted values, which could round the difference away. adjoint now uses exact integer checks or explicit bounds built from epsilon, dimension and operand magnitude. The 2×2 extreme-scale Bunch–Kaufman tolerance tightens from 1e-3 to 8 * 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_solverbase is skipped when rank()'s clamp at numeric_limits::min() drops a positive singular value \ge threshold() * sigma_0, because m \cdot \mathrm{solve}(m x) = m x cannot hold then. svd_least_square covers those inputs. Zero matrices keep the check (jacobisvd_56, bdcsvd_45).
  • cwise_min_max works around a GCC 10 AArch64 miscompile: the compiler kept the unreduced 128-bit partial result of maxCoeff() in v8 across a call, but AAPCS64 preserves only the low 64 bits.
  • New: approx_comparisons (four parts, covering FTZ, extreme scales, Homogeneous, lazy sums, complex-exponent pow, CustomComplex, and the no-allocation property), numerical_test_helpers, and benchmarks/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::isOrthogonal reports a nonzero vector as orthogonal to itself at scales 1e100 and 1e-100.
  • Sparse isApprox and complex scalar internal::isApprox accept twofold differences at scales 1e200 and 1e-200.
  • Real scalar internal::isApprox, which isApproxToConstant also 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.
Edited by Rasmus Munk Larsen

Merge request reports

Loading
Loading