SparseCore: Compare sparse matrices without squaring their norms libeigen/eigen!3353 Co-authored-by: Rasmus Munk Larsen <rmlarsen@gmail.com>
diff --git a/Eigen/src/SparseCore/BlockSparseMatrix.h b/Eigen/src/SparseCore/BlockSparseMatrix.h index f79d743..a9ef4b3 100644 --- a/Eigen/src/SparseCore/BlockSparseMatrix.h +++ b/Eigen/src/SparseCore/BlockSparseMatrix.h
@@ -601,16 +601,19 @@ // ------------------------------------------------------------------------- bool isApprox(const BlockSparseMatrix& other, const typename NumTraits<Scalar>::Real& prec = NumTraits<Scalar>::dummy_precision()) const { - using RealScalar = typename NumTraits<Scalar>::Real; + using Comparison = internal::sparse_approx_comparison<Scalar>; + using Wide = typename Comparison::Wide; // Frobenius-norm comparison, matching SparseMatrixBase::isApprox semantics // but computed directly from the block values — no scalar-level SparseMatrix - // materialization of either operand. Explicit zero blocks contribute 0 to - // every norm, so the result is independent of structural differences. - RealScalar n2a = m_values.head(nonZeros()).matrix().squaredNorm(); - RealScalar n2b = other.m_values.head(other.nonZeros()).matrix().squaredNorm(); - BlockSparseMatrix diff = *this - other; - RealScalar d2 = diff.m_values.head(diff.nonZeros()).matrix().squaredNorm(); - return d2 <= prec * prec * numext::mini(n2a, n2b); + // materialization of either operand while the squared norms are in range. + // Explicit zero blocks contribute 0 to every norm, so the result is + // independent of structural differences. + const auto squaredNorm = [](const BlockSparseMatrix& m) { + return m.m_values.head(m.nonZeros()).matrix().template cast<Wide>().squaredNorm(); + }; + return Comparison::compare( + squaredNorm(*this), squaredNorm(other), nonZeros() + other.nonZeros(), prec, + [&] { return squaredNorm(*this - other); }, [&] { return toSparse().isApprox(other.toSparse(), prec); }); } // -------------------------------------------------------------------------
diff --git a/Eigen/src/SparseCore/SparseApproxComparisons.h b/Eigen/src/SparseCore/SparseApproxComparisons.h index 8db78d6..2af41d2 100644 --- a/Eigen/src/SparseCore/SparseApproxComparisons.h +++ b/Eigen/src/SparseCore/SparseApproxComparisons.h
@@ -16,6 +16,123 @@ namespace Eigen { +namespace internal { + +template <typename Scalar, bool UseScaledComparison> +struct sparse_approx_comparison { + using RealScalar = typename NumTraits<Scalar>::Real; + using Accumulator = RealScalar; + // Custom scalars are summed as they are; a conversion to std::complex is not required of them. + using Wide = Scalar; + + template <typename DifferenceNorm2, typename Scaled> + static bool compare(const RealScalar& x2, const RealScalar& y2, Index, const RealScalar& prec, + const DifferenceNorm2& difference2, const Scaled&) { + return difference2() <= prec * prec * numext::mini(x2, y2); + } + + template <typename X, typename Y> + static bool isApprox(const X& x, const Y& y, const RealScalar& prec) { + return compare( + x.squaredNorm(), y.squaredNorm(), 0, prec, [&] { return (x - y).squaredNorm(); }, [] { return false; }); + } +}; + +// |x - y| <= prec min(|x|, |y|) as in the dense approx_comparison_impl: from the squared norms of the stored values +// while they are in range, and otherwise from norms of values scaled by powers of two, so tiny and huge matrices are +// compared rather than their flushed or overflowed squares. +template <typename Scalar> +struct sparse_approx_comparison<Scalar, true> { + using RealScalar = typename NumTraits<Scalar>::Real; + using Accumulator = typename stable_norm_accumulator<RealScalar>::type; + // Coefficients widened to the accumulator, as approx_comparison_wide_t widens an expression's. + using Wide = std::conditional_t<NumTraits<Scalar>::IsComplex, std::complex<Accumulator>, Accumulator>; + using Dense = approx_comparison_impl<Scalar, true>; + + // x2 and y2 are the squared norms of operands with at most size stored values between them; difference2() gives + // that of their difference, and scaled() decides the comparison when the squares cannot. + template <typename DifferenceNorm2, typename Scaled> + static bool compare(const Accumulator& x2, const Accumulator& y2, Index size, const RealScalar& prec, + const DifferenceNorm2& difference2, const Scaled& scaled) { + const Accumulator minimum = numext::mini(x2, y2); + const Accumulator precision2 = Accumulator(prec) * Accumulator(prec); + const Accumulator bound = precision2 * minimum; + // Only the smaller norm enters the bound; overflow of the larger norm is harmless. + if (Dense::safe_squared_norm(bound, size) && + minimum >= Dense::template squared_norm_lower_bound<RealScalar>(size) && + precision2 >= Dense::template squared_norm_lower_bound<RealScalar>(1)) + return Accumulator(difference2()) <= bound; + return scaled(); + } + + template <typename X, typename Y> + static bool isApprox(const X& x, const Y& y, const RealScalar& prec) { + // Widening must not admit operands that cannot be subtracted, such as half and float. + EIGEN_CHECK_BINARY_COMPATIBILITY(scalar_difference_op<typename X::Scalar EIGEN_COMMA typename Y::Scalar>, + typename X::Scalar, typename Y::Scalar) + const auto& wideX = x.template cast<approx_comparison_wide_t<X>>(); + const auto& wideY = y.template cast<approx_comparison_wide_t<Y>>(); + // x - y stores at most the entries of both operands. + return compare( + wideX.squaredNorm(), wideY.squaredNorm(), stored_size(x) + stored_size(y), prec, + [&] { return (wideX - wideY).squaredNorm(); }, [&] { return isApprox_scaled(x, y, prec); }); + } + + private: + // An upper bound on the number of stored values; only types with compressed access count them. + template <typename X, std::enable_if_t<has_compressed_access<X>::value, int> = 0> + static Index stored_size(const X& x) { + return x.nonZeros(); + } + template <typename X, std::enable_if_t<!has_compressed_access<X>::value, int> = 0> + static Index stored_size(const X& x) { + return x.rows() * x.cols(); + } + + // The norm of a compressed matrix is that of its stored values, which the dense scaled norm takes as a vector. + template <typename Plain> + static comparison_magnitude<Accumulator> norm(const Plain& m) { + return scaled_comparison_norm(m.coeffs().matrix()); + } + + template <typename X, typename Y> + static EIGEN_DONT_INLINE bool isApprox_scaled(const X& xExpr, const Y& yExpr, const RealScalar& prec) { + constexpr int Order = X::IsRowMajor ? RowMajor : ColMajor; + using StorageIndex = typename X::StorageIndex; + using WideX = approx_comparison_wide_t<X>; + using WideY = approx_comparison_wide_t<Y>; + using WideDifference = std::conditional_t<NumTraits<WideX>::IsComplex, WideX, WideY>; + const SparseMatrix<WideX, Order, StorageIndex> x = xExpr.template cast<WideX>(); + const SparseMatrix<WideY, Order, StorageIndex> y = yExpr.template cast<WideY>(); + const auto nx = norm(x); + const auto ny = norm(y); + if (!nx.isFinite() || !ny.isFinite()) return false; + auto tolerance = nx <= ny ? nx : ny; + tolerance.multiply(numext::abs(Accumulator(prec))); + + // As in scaled_comparison_distance_impl: lift a maximum below 1 by a power of two first, so FTZ does not flush a + // subnormal difference of normal values, and halve if the difference of finite values overflows. + const Accumulator maxCoeff = + numext::mini(Accumulator(1), numext::maxi(scaled_comparison_max_coeff(x.coeffs().matrix().realView()), + scaled_comparison_max_coeff(y.coeffs().matrix().realView()))); + const safe_scaling_factors<Accumulator> factors = supports_power_of_two_scaling<Accumulator>::value + ? safe_scaling<Accumulator>::compute_floor_factors(maxCoeff) + : safe_scaling_factors<Accumulator>(); + SparseMatrix<WideDifference, Order, StorageIndex> difference = x * factors.invScale - y * factors.invScale; + auto distance = norm(difference); + if (!distance.isFinite()) { + const Accumulator halfInvScale = factors.invScale * Accumulator(0.5); + difference = x * halfInvScale - y * halfInvScale; + distance = norm(difference); + distance.multiply(Accumulator(2)); + } + distance.multiply(factors.scale); + return distance <= tolerance; + } +}; + +} // namespace internal + template <typename Derived> template <typename OtherDerived> bool SparseMatrixBase<Derived>::isApprox(const SparseMatrixBase<OtherDerived>& other, const RealScalar& prec) const { @@ -24,7 +141,7 @@ const typename internal::nested_eval<OtherDerived, 2, PlainObject>::type, const PlainObject> actualB(other.derived()); - return (actualA - actualB).squaredNorm() <= prec * prec * numext::mini(actualA.squaredNorm(), actualB.squaredNorm()); + return internal::sparse_approx_comparison<Scalar>::isApprox(actualA, actualB, prec); } } // end namespace Eigen
diff --git a/Eigen/src/SparseCore/SparseUtil.h b/Eigen/src/SparseCore/SparseUtil.h index e44d197..f747028 100644 --- a/Eigen/src/SparseCore/SparseUtil.h +++ b/Eigen/src/SparseCore/SparseUtil.h
@@ -69,6 +69,10 @@ template <typename Derived> struct has_compressed_access : bool_constant<(int(traits<Derived>::Flags) & CompressedAccessBit) != 0> {}; +// Defined in SparseApproxComparisons.h; BlockSparseMatrix.h uses it before that. +template <typename Scalar, bool = use_scaled_comparison<Scalar>::value> +struct sparse_approx_comparison; + template <typename T, int Rows, int Cols, int Flags> struct sparse_eval;
diff --git a/test/block_sparse_matrix.cpp b/test/block_sparse_matrix.cpp index cb781d2..9bd76bb 100644 --- a/test/block_sparse_matrix.cpp +++ b/test/block_sparse_matrix.cpp
@@ -8,6 +8,7 @@ // SPDX-License-Identifier: MPL-2.0 #include "sparse.h" +#include "CustomComplex.h" using namespace Eigen; @@ -685,6 +686,53 @@ VERIFY_IS_EQUAL(A.coeff(1, 3), 5); } +// isApprox() compared squared norms, which overflow to inf <= inf or underflow to 0 <= 0, so A.isApprox(3 * A) held for +// tiny and huge matrices. +template <typename Scalar, int Options> +void test_block_sparse_isapprox_extreme_scales() { + using RealScalar = typename NumTraits<Scalar>::Real; + using BSM = BlockSparseMatrix<Scalar, Options, 2, 2>; + const RealScalar highest = NumTraits<RealScalar>::highest(); + const RealScalar min = (std::numeric_limits<RealScalar>::min)(); + const RealScalar scales[] = {RealScalar(1), + highest / RealScalar(2), + highest / RealScalar(8), + numext::sqrt(highest), + RealScalar(1) / numext::sqrt(highest), + min * RealScalar(16), + std::numeric_limits<RealScalar>::denorm_min() * RealScalar(64)}; + for (const RealScalar s : scales) { + SparseMatrix<Scalar, Options> sparse(4, 6); + sparse.insert(0, 0) = Scalar(s); + sparse.insert(1, 1) = Scalar(-2 * s); + sparse.insert(3, 4) = Scalar(s / 2); + sparse.makeCompressed(); + const BSM a = BSM::fromSparse(sparse); + VERIFY(a.isApprox(a)); + VERIFY(!a.isApprox(a * Scalar(3))); + VERIFY(!a.isApprox(-a)); + VERIFY(!a.isApprox(BSM(2, 3))); + } +} + +// Custom scalars keep the squared comparison, so isApprox() must not widen them to std::complex. +void test_block_sparse_isapprox_custom_complex() { + using Scalar = CustomComplex<double>; + using BSM = BlockSparseMatrix<Scalar, ColMajor, 2, 2>; + STATIC_CHECK(!internal::use_scaled_comparison<Scalar>::value); + SparseMatrix<Scalar> sparse(4, 6); + sparse.insert(0, 0) = Scalar(1, 2); + sparse.insert(1, 1) = Scalar(-2, 0.5); + sparse.insert(3, 4) = Scalar(0.5, -1); + sparse.makeCompressed(); + const BSM a = BSM::fromSparse(sparse); + VERIFY(a.isApprox(a)); + VERIFY(!a.isApprox(a * Scalar(3))); + VERIFY(!a.isApprox(-a)); + VERIFY(sparse.isApprox(sparse)); + VERIFY(!sparse.isApprox(-sparse)); +} + void test_block_triplet_traits() { // Flat scalar array means BlockTriplet should be trivially copyable and // standard-layout for any trivially-copyable Scalar and StorageIndex. @@ -765,6 +813,10 @@ // BlockTriplet type traits CALL_SUBTEST_17(test_block_triplet_traits()); + CALL_SUBTEST_17((test_block_sparse_isapprox_extreme_scales<float, ColMajor>())); + CALL_SUBTEST_17((test_block_sparse_isapprox_extreme_scales<double, RowMajor>())); + CALL_SUBTEST_17((test_block_sparse_isapprox_extreme_scales<std::complex<double>, ColMajor>())); + CALL_SUBTEST_17(test_block_sparse_isapprox_custom_complex()); // Complex scalar coverage: conjugation paths in adjoint, selfadjoint, triangular products CALL_SUBTEST_18((test_block_sparse<2, 2, ColMajor, std::complex<double>>(4, 6)));
diff --git a/test/sparse_basic.cpp b/test/sparse_basic.cpp index 5ce505d..8462185 100644 --- a/test/sparse_basic.cpp +++ b/test/sparse_basic.cpp
@@ -1334,6 +1334,59 @@ materialize_selfadjoint_pattern_random_impl<Upper>(); } +// isApprox() compared squared norms, which overflow to inf <= inf or underflow to 0 <= 0, so a.isApprox(3 * a) held for +// tiny and huge matrices, and a.isApprox(-a) where a - (-a) overflows. +template <typename Scalar, int Options> +void sparse_isapprox_extreme_scales() { + using RealScalar = typename NumTraits<Scalar>::Real; + using SpMat = SparseMatrix<Scalar, Options>; + using OtherSpMat = SparseMatrix<Scalar, Options ^ RowMajor>; + const RealScalar highest = NumTraits<RealScalar>::highest(); + const RealScalar min = (std::numeric_limits<RealScalar>::min)(); + // highest / 2 makes a - (-a) overflow although a and -a are finite. + const RealScalar scales[] = {RealScalar(1), + highest / RealScalar(2), + highest / RealScalar(8), + numext::sqrt(highest), + RealScalar(1) / numext::sqrt(highest), + min * RealScalar(16), + std::numeric_limits<RealScalar>::denorm_min() * RealScalar(64)}; + for (const RealScalar s : scales) { + SpMat a(5, 4); + a.insert(0, 0) = Scalar(s); + a.insert(2, 1) = Scalar(-2 * s); + a.insert(4, 3) = Scalar(s / 2); + a.makeCompressed(); + const SpMat b = Scalar(3) * a; + SpMat uncompressed = a; + uncompressed.uncompress(); + const OtherSpMat other = a; + + VERIFY(a.isApprox(a)); + VERIFY(!a.isApprox(b)); + VERIFY(!b.isApprox(a)); + VERIFY(!a.isApprox(SpMat(-a))); + VERIFY(!a.isApprox(SpMat(a.rows(), a.cols()))); + VERIFY(uncompressed.isApprox(a)); + VERIFY(!uncompressed.isApprox(b)); + VERIFY(other.isApprox(a)); + VERIFY(!a.isApprox(Scalar(3) * other)); + + SparseVector<Scalar> v(6); + v.insert(1) = Scalar(s); + v.insert(4) = Scalar(-s); + VERIFY(v.isApprox(v)); + VERIFY(!v.isApprox(SparseVector<Scalar>(Scalar(3) * v))); + + // A tight prec: |a - c| = 2^-10 |a| exactly where the values stay normal. + if (s >= min * RealScalar(16) && s <= highest / RealScalar(8)) { + const SpMat c = a * Scalar(RealScalar(1) + RealScalar(1) / RealScalar(1024)); + VERIFY(a.isApprox(c, RealScalar(1) / RealScalar(512))); + VERIFY(!a.isApprox(c, RealScalar(1) / RealScalar(2048))); + } + } +} + template <int> void bug1105() { // Regression test for bug 1105 @@ -1383,6 +1436,9 @@ CALL_SUBTEST_5((big_sparse_triplet<SparseMatrix<double, ColMajor, long int>>(10000, 10000, 0.125))); CALL_SUBTEST_5(bug1105<0>()); + CALL_SUBTEST_1((sparse_isapprox_extreme_scales<float, ColMajor>())); + CALL_SUBTEST_1((sparse_isapprox_extreme_scales<double, RowMajor>())); + CALL_SUBTEST_1((sparse_isapprox_extreme_scales<std::complex<double>, ColMajor>())); CALL_SUBTEST_1(sparse_sub_assign_eigenbase<0>()); CALL_SUBTEST_1(ambivector_coeff<0>()); CALL_SUBTEST_1(sparsity_pattern_ref_sparse_expressions<0>());