diff --git a/cpp/tests/sparse/solver/lanczos.cu b/cpp/tests/sparse/solver/lanczos.cu index cd10cee71a..c66aa2f8e3 100644 --- a/cpp/tests/sparse/solver/lanczos.cu +++ b/cpp/tests/sparse/solver/lanczos.cu @@ -8,9 +8,13 @@ #include #include #include +#include #include #include +#include #include +#include +#include #include #include #include @@ -35,7 +39,9 @@ #include #include +#include #include +#include namespace raft::sparse::solver { @@ -52,7 +58,6 @@ struct lanczos_inputs { std::vector rows; // indptr std::vector cols; // indices std::vector vals; // data - std::vector expected_eigenvalues; }; template @@ -68,9 +73,175 @@ struct rmat_lanczos_inputs { int r_scale; int c_scale; float sparsity; - std::vector expected_eigenvalues; }; +/** + * @brief Compute the Frobenius norm of a sparse matrix from its stored values. + * + * Used to scale residual/orthogonality tolerances to the magnitude of the + * problem instead of an absolute value. + * + * @tparam ValueType data type of the matrix values (float or double) + * @tparam NnzType integral type of the nonzero count + * @param[in] handle raft resources + * @param[in] vals device pointer to the nnz stored values of the matrix + * @param[in] nnz number of stored values + * @return the Frobenius norm sqrt(sum(vals[i]^2)) + */ +template +ValueType compute_frobenius_norm(raft::resources const& handle, ValueType const* vals, NnzType nnz) +{ + auto vals_view = + raft::make_device_vector_view(vals, static_cast(nnz)); + ValueType sum_sq{}; + raft::linalg::dot(handle, vals_view, vals_view, raft::make_host_scalar_view(&sum_sq)); + raft::resource::sync_stream(handle); + return std::sqrt(sum_sq); +} + +/** + * @brief Validate that (eigenvalues[i], eigenvectors[:, i]) are genuine + * eigenpairs of `A`, without relying on hardcoded golden eigenvalues. + * + * Checks the defining characteristics of a converged real-symmetric + * eigensolve: + * - ||A*v_i - lambda_i*v_i|| is within the noise floor implied by the + * solver's own convergence criterion (see lanczos_solver_config::tolerance, + * which is documented as governing exactly this residual). + * - lambda_i matches the Rayleigh quotient of v_i. This pins the eigenvalue + * down more sharply than the residual alone on matrices with a large + * ||A||, where a small eigenvalue's residual budget is dominated by the + * matrix scale. + * - each eigenvector is unit-normalized. + * - distinct eigenvectors are mutually orthogonal (catches "ghost" + * duplicate eigenpairs from loss of orthogonality, which a residual-only + * check would miss). + * - eigenvalues are returned in ascending algebraic order (verified + * against every previously-hardcoded golden array in this file, which + * were all ascending regardless of `which`). + * + * @tparam IndexType integral type of the matrix indices + * @tparam ValueType data type of matrix values and eigenpairs (float or double) + * @param[in] handle raft resources + * @param[in] A sparse matrix wrapper providing the A*v product + * @param[in] eigenvalues device pointer to n_components computed eigenvalues + * @param[in] eigenvectors device pointer to the column-major n x n_components + * eigenvector matrix + * @param[in] n dimension of the (square) matrix + * @param[in] n_components number of eigenpairs to validate + * @param[in] frobenius_norm Frobenius norm of A, used to scale tolerances + */ +template +void expect_valid_eigenpairs(raft::resources const& handle, + raft::spectral::matrix::sparse_matrix_t const& A, + ValueType const* eigenvalues, + ValueType* eigenvectors, + IndexType n, + int n_components, + ValueType frobenius_norm) +{ + auto stream = resource::get_cuda_stream(handle); + auto eps = std::numeric_limits::epsilon(); + + // Tolerances are expressed as multiples of the floating-point noise floor + // for each kind of quantity, since the exact low-order bits depend on + // cuSPARSE/cuBLAS reduction order, which varies by CUDA version and GPU + // architecture (the root cause of issue #2519). + // + // Quantities carrying the magnitude of A (the residual ||A*v - lambda*v|| + // and the Rayleigh quotient) have a noise floor proportional to + // ||A||_F * eps. Measured across all fixtures in this file on sm_75 / + // CUDA 13, these stay within 15 * ||A||_F * eps and show no dependence on + // n (n=100 and n=4096 both land in that range), so n is deliberately not a + // factor here -- including it would inflate the tolerance ~40x on the RMAT + // fixture and let a 1% eigenvalue error pass undetected. + // + // Dimensionless quantities (deviation from unit norm, mutual + // orthogonality) have a noise floor of a small multiple of eps; measured + // values stay within 100 * eps. + const ValueType matrix_tol = ValueType(500) * std::max(frobenius_norm, ValueType(1)) * eps; + const ValueType unit_tol = ValueType(1000) * eps; + const ValueType norm_tol = unit_tol; + const ValueType ortho_tol = unit_tol; + const ValueType residual_tol = matrix_tol; + + std::vector host_eigenvalues(n_components); + raft::update_host(host_eigenvalues.data(), eigenvalues, n_components, stream); + raft::resource::sync_stream(handle); + + for (int i = 1; i < n_components; ++i) { + EXPECT_LE(host_eigenvalues[i - 1], host_eigenvalues[i] + matrix_tol) + << "eigenvalues not returned in ascending algebraic order at index " << i; + } + + auto residual = raft::make_device_vector(handle, n); + + for (int i = 0; i < n_components; ++i) { + ValueType* v_i = eigenvectors + static_cast(i) * n; + auto v_i_view = raft::make_device_vector_view(v_i, n); + + // residual = A * v_i + A.mv(ValueType{1}, v_i, ValueType{0}, residual.data_handle()); + + ValueType rayleigh{}; + raft::linalg::dot(handle, + v_i_view, + raft::make_const_mdspan(residual.view()), + raft::make_host_scalar_view(&rayleigh)); + + // residual = A*v_i - lambda_i*v_i + ValueType neg_lambda = -host_eigenvalues[i]; + raft::linalg::axpy( + handle, raft::make_host_scalar_view(&neg_lambda), v_i_view, residual.view()); + + ValueType residual_sq{}; + raft::linalg::dot(handle, + raft::make_const_mdspan(residual.view()), + raft::make_const_mdspan(residual.view()), + raft::make_host_scalar_view(&residual_sq)); + + ValueType v_norm_sq{}; + raft::linalg::dot( + handle, v_i_view, v_i_view, raft::make_host_scalar_view(&v_norm_sq)); + + // dot with a host-scalar output is not guaranteed to have synchronized; + // make the host-side reads below safe. + raft::resource::sync_stream(handle); + + ValueType residual_norm = std::sqrt(residual_sq); + ValueType v_norm = std::sqrt(v_norm_sq); + + EXPECT_NEAR(v_norm, ValueType(1), norm_tol) + << "eigenvector " << i << " is not unit-normalized: ||v||=" << v_norm; + EXPECT_LE(residual_norm, residual_tol) + << "eigenpair " << i + << " failed ||A*v - lambda*v|| residual check: residual=" << residual_norm + << " lambda=" << host_eigenvalues[i] << " tol=" << residual_tol; + + // The Rayleigh quotient (v^T A v) / (v^T v) is the eigenvalue implied by + // the returned eigenvector. Comparing it against the returned eigenvalue + // pins down lambda independently of the residual check, which on a + // large-||A|| matrix is too coarse to notice a small relative error in a + // small-magnitude eigenvalue. + EXPECT_NEAR(host_eigenvalues[i], rayleigh / v_norm_sq, matrix_tol) + << "eigenvalue " << i + << " disagrees with the Rayleigh quotient of its eigenvector: lambda=" << host_eigenvalues[i] + << " rayleigh=" << rayleigh / v_norm_sq; + + for (int j = 0; j < i; ++j) { + ValueType const* v_j = eigenvectors + static_cast(j) * n; + auto v_j_view = raft::make_device_vector_view(v_j, n); + ValueType dot_ij{}; + raft::linalg::dot( + handle, v_i_view, v_j_view, raft::make_host_scalar_view(&dot_ij)); + raft::resource::sync_stream(handle); + EXPECT_LE(std::abs(dot_ij), ortho_tol) + << "eigenvectors " << i << " and " << j + << " are not orthogonal: |v_i . v_j|=" << std::abs(dot_ij); + } + } +} + template class rmat_lanczos_tests : public ::testing::TestWithParam> { @@ -79,8 +250,6 @@ class rmat_lanczos_tests : params(::testing::TestWithParam>::GetParam()), stream(resource::get_cuda_stream(handle)), rng(params.seed), - expected_eigenvalues(raft::make_device_vector( - handle, params.n_components)), r_scale(params.r_scale), c_scale(params.c_scale), sparsity(params.sparsity) @@ -88,13 +257,7 @@ class rmat_lanczos_tests } protected: - void SetUp() override - { - raft::copy(expected_eigenvalues.data_handle(), - params.expected_eigenvalues.data(), - params.n_components, - stream); - } + void SetUp() override {} void TearDown() override {} @@ -200,11 +363,15 @@ class rmat_lanczos_tests eigenvalues.view(), eigenvectors.view()); - ASSERT_TRUE(raft::devArrMatch(eigenvalues.data_handle(), - expected_eigenvalues.data_handle(), - n_components, - raft::CompareApprox(1e-5), - stream)); + ValueType frobenius_norm = + compute_frobenius_norm(handle, symmetric_coo.vals(), symmetric_coo.nnz); + expect_valid_eigenpairs(handle, + csr_m, + eigenvalues.data_handle(), + eigenvectors.data_handle(), + (IndexType)symmetric_coo.n_rows, + n_components, + frobenius_norm); // Reproducibility test - run again with same seed and verify exact match raft::device_vector eigenvalues2 = @@ -254,11 +421,13 @@ class rmat_lanczos_tests eigenvalues_coo.view(), eigenvectors_coo.view()); - ASSERT_TRUE(raft::devArrMatch(eigenvalues_coo.data_handle(), - expected_eigenvalues.data_handle(), - n_components, - raft::CompareApprox(1e-4), - stream)); + expect_valid_eigenpairs(handle, + csr_m, + eigenvalues_coo.data_handle(), + eigenvectors_coo.data_handle(), + (IndexType)symmetric_coo.n_rows, + n_components, + frobenius_norm); } protected: @@ -269,7 +438,6 @@ class rmat_lanczos_tests int r_scale; int c_scale; float sparsity; - raft::device_vector expected_eigenvalues; }; template @@ -288,9 +456,7 @@ class lanczos_tests : public ::testing::TestWithParam( handle, params.n_components)), eigenvectors(raft::make_device_matrix( - handle, n, params.n_components)), - expected_eigenvalues( - raft::make_device_vector(handle, params.n_components)) + handle, n, params.n_components)) { } @@ -300,10 +466,6 @@ class lanczos_tests : public ::testing::TestWithParam( const_cast(vals.data_handle()), csr_structure); + raft::spectral::matrix::sparse_matrix_t const A_check{ + handle, rows.data_handle(), cols.data_handle(), vals.data_handle(), n, (uint64_t)nnz}; + ValueType frobenius_norm = compute_frobenius_norm(handle, vals.data_handle(), nnz); + std::get<0>(stats) = raft::sparse::solver::lanczos_compute_eigenpairs( handle, config, @@ -348,13 +514,13 @@ class lanczos_tests : public ::testing::TestWithParam( - eigenvalues.data_handle(), - expected_eigenvalues.data_handle(), - params.n_components, - raft::CompareApprox( - params.which == raft::sparse::solver::LANCZOS_WHICH::SM ? 5e-5 : 1e-5), - stream)); + expect_valid_eigenpairs(handle, + A_check, + eigenvalues.data_handle(), + eigenvectors.data_handle(), + (IndexType)n, + params.n_components, + frobenius_norm); // Reproducibility test - run again with same seed and verify exact match raft::device_vector eigenvalues2 = @@ -403,13 +569,13 @@ class lanczos_tests : public ::testing::TestWithParam( - eigenvalues_coo.data_handle(), - expected_eigenvalues.data_handle(), - params.n_components, - raft::CompareApprox( - params.which == raft::sparse::solver::LANCZOS_WHICH::SM ? 5e-5 : 1e-5), - stream)); + expect_valid_eigenpairs(handle, + A_check, + eigenvalues_coo.data_handle(), + eigenvectors_coo.data_handle(), + (IndexType)n, + params.n_components, + frobenius_norm); } protected: @@ -425,7 +591,6 @@ class lanczos_tests : public ::testing::TestWithParam v0; raft::device_vector eigenvalues; raft::device_matrix eigenvectors; - raft::device_vector expected_eigenvalues; }; template @@ -695,8 +860,7 @@ const std::vector> inputsf = { 0.4350971, 0.6997072, 0.4320931, 0.3315690, 0.0844443, 0.1445242, 0.3059566, 0.6594226, 0.8961608, 0.6498466, 0.9585592, 0.7827352, 0.6498466, 0.2812338, 0.1767728, 0.5810611, 0.7269946, 0.6997072, 0.1705930, 0.1792683, 0.1077409, 0.9368132, 0.4823034, 0.8311127, - 0.7194629, 0.6273088, 0.2909178, 0.5188584, 0.5876446, 0.2812338}, - {-2.0369630, -1.7673520}}}; + 0.7194629, 0.6273088, 0.2909178, 0.5188584, 0.5876446, 0.2812338}}}; const std::vector> inputsd = { {2, @@ -746,8 +910,7 @@ const std::vector> inputsd = { 0.4350971, 0.6997072, 0.4320931, 0.3315690, 0.0844443, 0.1445242, 0.3059566, 0.6594226, 0.8961608, 0.6498466, 0.9585592, 0.7827352, 0.6498466, 0.2812338, 0.1767728, 0.5810611, 0.7269946, 0.6997072, 0.1705930, 0.1792683, 0.1077409, 0.9368132, 0.4823034, 0.8311127, - 0.7194629, 0.6273088, 0.2909178, 0.5188584, 0.5876446, 0.2812338}, - {-2.0369630, -1.7673520}}}; + 0.7194629, 0.6273088, 0.2909178, 0.5188584, 0.5876446, 0.2812338}}}; const std::vector> inputsd_SM = { {2, @@ -760,8 +923,7 @@ const std::vector> inputsd_SM = { 42, rows(), cols(), - vals(), - {-0.03944135, 0.01367824}}}; + vals()}}; const std::vector> inputsd_LM = { {2, @@ -774,8 +936,7 @@ const std::vector> inputsd_LM = { 42, rows(), cols(), - vals(), - {-2.00968758, 3.45939575}}}; + vals()}}; const std::vector> inputsd_LA = { {2, @@ -788,8 +949,7 @@ const std::vector> inputsd_LA = { 42, rows(), cols(), - vals(), - {1.99482678, 3.45939575}}}; + vals()}}; const std::vector> inputsd_SA = { {2, @@ -802,8 +962,7 @@ const std::vector> inputsd_SA = { 42, rows(), cols(), - vals(), - {-2.00968758, -1.83487935}}}; + vals()}}; const std::vector> inputsf_SM = { {2, @@ -816,8 +975,7 @@ const std::vector> inputsf_SM = { 42, rows(), cols(), - vals(), - {-0.03944135, 0.01367824}}}; + vals()}}; const std::vector> inputsf_LM = { {2, @@ -830,8 +988,7 @@ const std::vector> inputsf_LM = { 42, rows(), cols(), - vals(), - {-2.00968758, 3.45939575}}}; + vals()}}; const std::vector> inputsf_LA = { {2, @@ -844,8 +1001,7 @@ const std::vector> inputsf_LA = { 42, rows(), cols(), - vals(), - {1.99482678, 3.45939575}}}; + vals()}}; const std::vector> inputsf_SA = { {2, @@ -858,29 +1014,10 @@ const std::vector> inputsf_SA = { 42, rows(), cols(), - vals(), - {-2.00968758, -1.83487935}}}; + vals()}}; const std::vector> rmat_inputsf = { - {50, - 100, - 10000, - 0, - 0, - 1e-9, - raft::sparse::solver::LANCZOS_WHICH::SA, - 42, - 12, - 12, - 1, - {-122.526794, -74.00686, -59.698284, -54.68617, -49.686813, -34.02644, -32.130703, - -31.26906, -30.32097, -22.946098, -20.497862, -20.23817, -19.269697, -18.42496, - -17.675667, -17.013401, -16.734581, -15.820215, -15.73925, -15.448187, -15.044634, - -14.692028, -14.127425, -13.967386, -13.6237755, -13.469393, -13.181225, -12.777589, - -12.623185, -12.55508, -12.2874565, -12.053391, -11.677346, -11.558279, -11.163732, - -10.922034, -10.7936945, -10.558049, -10.205776, -10.005316, -9.559181, -9.491834, - -9.242631, -8.883637, -8.765364, -8.688508, -8.458255, -8.385196, -8.217982, - -8.0442095}}}; + {50, 100, 10000, 0, 0, 1e-9, raft::sparse::solver::LANCZOS_WHICH::SA, 42, 12, 12, 1}}; using LanczosTestF = lanczos_tests; TEST_P(LanczosTestF, Result) { Run(); }