diff --git a/cpp/include/raft/sparse/linalg/detail/laplacian.cuh b/cpp/include/raft/sparse/linalg/detail/laplacian.cuh index 3325c9717d..fc6ef79ad3 100644 --- a/cpp/include/raft/sparse/linalg/detail/laplacian.cuh +++ b/cpp/include/raft/sparse/linalg/detail/laplacian.cuh @@ -1,13 +1,20 @@ /* - * SPDX-FileCopyrightText: Copyright (c) 2025, NVIDIA CORPORATION. + * SPDX-FileCopyrightText: Copyright (c) 2025-2026, NVIDIA CORPORATION. * SPDX-License-Identifier: Apache-2.0 */ #pragma once #include #include #include +#include +#include #include #include +#include + +#include +#include +#include #include @@ -108,6 +115,123 @@ auto compute_graph_laplacian( return result; } +template +device_coo_matrix compute_graph_laplacian( + raft::resources const& res, device_coo_matrix_view input) +{ + auto input_structure = input.structure_view(); + auto dim = input_structure.get_n_rows(); + RAFT_EXPECTS(dim == input_structure.get_n_cols(), + "The graph Laplacian can only be computed on a square adjacency matrix"); + + auto stream = resource::get_cuda_stream(res); + + auto marked_diagonal = raft::make_device_vector(res, dim); + raft::matrix::fill(res, marked_diagonal.view(), int(1)); + auto marked_diagonal_ptr = marked_diagonal.data_handle(); + auto rows_ptr = input_structure.get_rows().data(); + auto cols_ptr = input_structure.get_cols().data(); + auto values_ptr = input.get_elements().data(); + + auto diagonal_count = raft::make_device_scalar(res, RowType(0)); + auto diagonal_count_ptr = diagonal_count.data_handle(); + + // mark which diagonal elements are present in the input matrix and also keep track of counter + raft::linalg::map_offset( + res, + raft::make_device_vector_view(values_ptr, input_structure.get_nnz()), + [rows_ptr, cols_ptr, values_ptr, marked_diagonal_ptr, diagonal_count_ptr] __device__(auto idx) { + if (rows_ptr[idx] == cols_ptr[idx]) { + marked_diagonal_ptr[rows_ptr[idx]] = 0; + atomicAdd(diagonal_count_ptr, RowType(1)); + } + return values_ptr[idx]; + }); + + auto host_diagonal_count = raft::make_host_scalar(RowType(0)); + raft::copy(host_diagonal_count.data_handle(), diagonal_count_ptr, 1, stream); + resource::sync_stream(res); + RowType extra_diagonal_space = dim - host_diagonal_count(0); + + auto result = make_device_coo_matrix, + std::remove_const_t, + std::remove_const_t, + std::remove_const_t>( + res, dim, dim, input_structure.get_nnz() + extra_diagonal_space); + + auto result_rows_ptr = result.structure_view().get_rows().data(); + auto result_cols_ptr = result.structure_view().get_cols().data(); + auto result_values_ptr = result.get_elements().data(); + + raft::copy(result_values_ptr, values_ptr, input_structure.get_nnz(), stream); + raft::copy(result_cols_ptr, cols_ptr, input_structure.get_nnz(), stream); + raft::copy(result_rows_ptr, rows_ptr, input_structure.get_nnz(), stream); + + auto scan_diagonal = raft::make_device_vector(res, dim); + auto scan_diagonal_ptr = scan_diagonal.data_handle(); + auto input_nnz = input_structure.get_nnz(); + + thrust::exclusive_scan(raft::resource::get_thrust_policy(res), + marked_diagonal_ptr, + marked_diagonal_ptr + dim, + scan_diagonal_ptr); + + // populate the extra diagonal indexes and initialize the values to 0 + raft::linalg::map_offset(res, + marked_diagonal.view(), + [result_rows_ptr, + result_cols_ptr, + result_values_ptr, + marked_diagonal_ptr, + scan_diagonal_ptr, + input_nnz] __device__(auto idx) { + if (marked_diagonal_ptr[idx] == 1) { + result_rows_ptr[input_nnz + scan_diagonal_ptr[idx]] = idx; + result_cols_ptr[input_nnz + scan_diagonal_ptr[idx]] = idx; + result_values_ptr[input_nnz + scan_diagonal_ptr[idx]] = 0; + } + return marked_diagonal_ptr[idx]; + }); + + raft::sparse::op::coo_sort( + dim, + dim, + result.structure_view().get_nnz(), + result.structure_view().get_rows().data(), + result.structure_view().get_cols().data(), + result.get_elements().data(), + raft::resource::get_cuda_stream(res)); + + auto result_nnz = result.structure_view().get_nnz(); + auto degrees = raft::make_device_vector(res, dim); + raft::matrix::fill(res, degrees.view(), ElementType(0)); + auto degrees_ptr = degrees.data_handle(); + + // D + thrust::reduce_by_key(raft::resource::get_thrust_policy(res), + result_rows_ptr, // keys_first + result_rows_ptr + result_nnz, // keys_last + result_values_ptr, // values_first + thrust::make_discard_iterator(), // keys_output (discarded) + degrees_ptr); // values_output (row sums) + + // D - A + raft::linalg::map_offset( + res, + raft::make_device_vector_view(result_values_ptr, result_nnz), + [result_rows_ptr, result_cols_ptr, result_values_ptr, degrees_ptr] __device__(auto idx) { + if (result_rows_ptr[idx] == result_cols_ptr[idx]) { + result_values_ptr[idx] = + degrees_ptr[result_rows_ptr[idx]] - result_values_ptr[idx]; // on diagonal + } else { + result_values_ptr[idx] = -result_values_ptr[idx]; // off diagonal + } + return result_values_ptr[idx]; + }); + + return result; +} + /** * @brief Given a CSR adjacency matrix, return the normalized graph Laplacian * @@ -129,11 +253,10 @@ auto compute_graph_laplacian( * * @return A CSR matrix containing the normalized graph Laplacian */ -template -auto laplacian_normalized( - raft::resources const& res, - device_csr_matrix_view input, - device_vector_view diagonal_out) +template +auto laplacian_normalized(raft::resources const& res, + SparseInputType input, + device_vector_view diagonal_out) { auto laplacian = detail::compute_graph_laplacian(res, input); auto laplacian_structure = laplacian.structure_view(); diff --git a/cpp/include/raft/sparse/linalg/laplacian.cuh b/cpp/include/raft/sparse/linalg/laplacian.cuh index 6f807d5b89..5b43b2e4d5 100644 --- a/cpp/include/raft/sparse/linalg/laplacian.cuh +++ b/cpp/include/raft/sparse/linalg/laplacian.cuh @@ -1,5 +1,5 @@ /* - * SPDX-FileCopyrightText: Copyright (c) 2025, NVIDIA CORPORATION. + * SPDX-FileCopyrightText: Copyright (c) 2025-2026, NVIDIA CORPORATION. * SPDX-License-Identifier: Apache-2.0 */ #pragma once @@ -22,6 +22,17 @@ auto compute_graph_laplacian( return detail::compute_graph_laplacian(res, input); } +/** Given a COO adjacency matrix, return the graph Laplacian + * + * Note that for non-symmetric matrices, the out-degree Laplacian is returned. + */ +template +device_coo_matrix compute_graph_laplacian( + raft::resources const& res, device_coo_matrix_view input) +{ + return detail::compute_graph_laplacian(res, input); +} + /** * @brief Given a CSR adjacency matrix, return the normalized graph Laplacian * @@ -49,7 +60,41 @@ auto laplacian_normalized( device_csr_matrix_view input, device_vector_view diagonal_out) { - return detail::laplacian_normalized(res, input, diagonal_out); + return detail::laplacian_normalized< + device_csr_matrix_view, + ElementType, + IndptrType>(res, input, diagonal_out); +} + +/** + * @brief Given a COO adjacency matrix, return the normalized graph Laplacian + * + * Return the normalized Laplacian matrix, which is defined as D^(-1/2) * L * D^(-1/2), + * where D is the diagonal degree matrix and L is the graph Laplacian. + * Also returns the scaled diagonal degree matrix. + * + * + * @tparam ElementType The data type of the matrix elements + * @tparam RowType The data type of the row pointers + * @tparam ColType The data type of the column indices + * @tparam NZType The data type for representing nonzero counts + * + * @param[in] res RAFT resources for managing device memory and streams + * @param[in] input View of the input COO adjacency matrix + * @param[out] diagonal_out View of the output vector where the scaled diagonal degree + * matrix D^(-1/2) will be stored (must be pre-allocated with + * size at least n_rows) + * + * @return A COO matrix containing the normalized graph Laplacian + */ +template +auto laplacian_normalized(raft::resources const& res, + device_coo_matrix_view input, + device_vector_view diagonal_out) +{ + return detail::laplacian_normalized, + ElementType, + RowType>(res, input, diagonal_out); } } // namespace raft::sparse::linalg diff --git a/cpp/include/raft/sparse/solver/lanczos.cuh b/cpp/include/raft/sparse/solver/lanczos.cuh index 3a5e244923..902396298e 100644 --- a/cpp/include/raft/sparse/solver/lanczos.cuh +++ b/cpp/include/raft/sparse/solver/lanczos.cuh @@ -1,5 +1,5 @@ /* - * SPDX-FileCopyrightText: Copyright (c) 2022-2025, NVIDIA CORPORATION. + * SPDX-FileCopyrightText: Copyright (c) 2022-2026, NVIDIA CORPORATION. * SPDX-License-Identifier: Apache-2.0 */ #ifndef __LANCZOS_H @@ -29,11 +29,11 @@ namespace raft::sparse::solver { * @param eigenvectors output eigenvectors * @return Zero if successful. Otherwise non-zero. */ -template +template auto lanczos_compute_smallest_eigenvectors( raft::resources const& handle, lanczos_solver_config const& config, - raft::device_csr_matrix_view A, + raft::device_csr_matrix_view A, std::optional> v0, raft::device_vector_view eigenvalues, raft::device_matrix_view eigenvectors) -> int @@ -54,11 +54,11 @@ auto lanczos_compute_smallest_eigenvectors( * @param eigenvectors output eigenvectors * @return Zero if successful. Otherwise non-zero. */ -template +template auto lanczos_compute_smallest_eigenvectors( raft::resources const& handle, lanczos_solver_config const& config, - raft::device_coo_matrix_view A, + raft::device_coo_matrix_view A, std::optional> v0, raft::device_vector_view eigenvalues, raft::device_matrix_view eigenvectors) -> int diff --git a/cpp/tests/sparse/laplacian.cu b/cpp/tests/sparse/laplacian.cu index 7825b3178b..03193cfca8 100644 --- a/cpp/tests/sparse/laplacian.cu +++ b/cpp/tests/sparse/laplacian.cu @@ -1,12 +1,17 @@ /* - * SPDX-FileCopyrightText: Copyright (c) 2020-2024, NVIDIA CORPORATION. + * SPDX-FileCopyrightText: Copyright (c) 2020-2026, NVIDIA CORPORATION. * SPDX-License-Identifier: Apache-2.0 */ +#include #include #include #include +#include +#include +#include #include +#include #include #include @@ -125,7 +130,7 @@ INSTANTIATE_TEST_SUITE_P(LaplacianTests, return info.param.name; }); -TEST(Raft, ComputeGraphLaplacianNormalized) +TEST(Raft, ComputeGraphLaplacianNormalizedCSR) { // Using the same adjacency matrix as in the ComputeGraphLaplacian test: // [[0 1 1 1] @@ -233,4 +238,191 @@ TEST(Raft, ComputeGraphLaplacianNormalized) EXPECT_EQ(expected_indptr, normalized_laplacian_indptr); } +TEST(Raft, ComputeGraphLaplacianNormalizedCOO) +{ + // Using the same adjacency matrix as in the ComputeGraphLaplacian test: + // [[0 1 1 1] + // [1 0 0 1] + // [1 0 0 0] + // [1 1 0 0]] + + auto data = std::vector{1, 1, 1, 1, 1, 1, 1, 1}; + auto indices = std::vector{1, 2, 3, 0, 3, 0, 0, 1}; + auto indptr = std::vector{0, 3, 5, 6, 8}; + + auto res = raft::resources{}; + + // First create a CSR matrix + auto adjacency_matrix_csr = make_device_csr_matrix( + res, (indptr.size() - 1), (indptr.size() - 1), data.size()); + auto adjacency_structure_csr = adjacency_matrix_csr.structure_view(); + raft::copy(adjacency_matrix_csr.get_elements().data(), + &(data[0]), + data.size(), + raft::resource::get_cuda_stream(res)); + raft::copy(adjacency_structure_csr.get_indices().data(), + &(indices[0]), + indices.size(), + raft::resource::get_cuda_stream(res)); + raft::copy(adjacency_structure_csr.get_indptr().data(), + &(indptr[0]), + indptr.size(), + raft::resource::get_cuda_stream(res)); + + // Convert CSR to COO + auto adjacency_matrix_coo = + make_device_coo_matrix(res, + adjacency_structure_csr.get_n_rows(), + adjacency_structure_csr.get_n_cols(), + adjacency_structure_csr.get_nnz()); + + raft::sparse::convert::csr_to_coo(adjacency_structure_csr.get_indptr().data(), + adjacency_structure_csr.get_n_rows(), + adjacency_matrix_coo.structure_view().get_rows().data(), + adjacency_structure_csr.get_nnz(), + raft::resource::get_cuda_stream(res)); + + raft::copy(adjacency_matrix_coo.structure_view().get_cols().data(), + adjacency_structure_csr.get_indices().data(), + adjacency_structure_csr.get_nnz(), + raft::resource::get_cuda_stream(res)); + + raft::copy(adjacency_matrix_coo.get_elements().data(), + adjacency_matrix_csr.get_elements().data(), + adjacency_structure_csr.get_nnz(), + raft::resource::get_cuda_stream(res)); + + // Create diagonal output vector + auto diagonal_out = + raft::make_device_vector(res, adjacency_structure_csr.get_n_rows()); + + // Compute normalized Laplacian using COO matrix (result is also COO) + auto normalized_laplacian_coo = + laplacian_normalized(res, adjacency_matrix_coo.view(), diagonal_out.view()); + auto normalized_laplacian_coo_structure = normalized_laplacian_coo.structure_view(); + + // Sort the COO matrix first + raft::sparse::op::coo_sort(normalized_laplacian_coo_structure.get_n_rows(), + normalized_laplacian_coo_structure.get_n_cols(), + normalized_laplacian_coo_structure.get_nnz(), + normalized_laplacian_coo_structure.get_rows().data(), + normalized_laplacian_coo_structure.get_cols().data(), + normalized_laplacian_coo.get_elements().data(), + raft::resource::get_cuda_stream(res)); + + // Convert COO result to CSR for comparison + auto normalized_laplacian_csr = + make_device_csr_matrix(res, + normalized_laplacian_coo_structure.get_n_rows(), + normalized_laplacian_coo_structure.get_n_cols(), + normalized_laplacian_coo_structure.get_nnz()); + auto normalized_laplacian_csr_structure = normalized_laplacian_csr.structure_view(); + + // Initialize indptr to zeros + raft::matrix::fill( + res, + raft::make_device_vector_view(normalized_laplacian_csr_structure.get_indptr().data(), + normalized_laplacian_coo_structure.get_n_rows() + 1), + int(0)); + + // Convert sorted COO to CSR + raft::sparse::convert::sorted_coo_to_csr( + normalized_laplacian_coo_structure.get_rows().data(), + normalized_laplacian_coo_structure.get_nnz(), + normalized_laplacian_csr_structure.get_indptr().data(), + normalized_laplacian_coo_structure.get_n_rows(), + raft::resource::get_cuda_stream(res)); + + // Manually set the last element of indptr to nnz (workaround for potential bug) + int nnz = normalized_laplacian_coo_structure.get_nnz(); + raft::copy(normalized_laplacian_csr_structure.get_indptr().data() + + normalized_laplacian_coo_structure.get_n_rows(), + &nnz, + 1, + raft::resource::get_cuda_stream(res)); + + raft::copy(normalized_laplacian_csr_structure.get_indices().data(), + normalized_laplacian_coo_structure.get_cols().data(), + normalized_laplacian_coo_structure.get_nnz(), + raft::resource::get_cuda_stream(res)); + + raft::copy(normalized_laplacian_csr.get_elements().data(), + normalized_laplacian_coo.get_elements().data(), + normalized_laplacian_coo_structure.get_nnz(), + raft::resource::get_cuda_stream(res)); + + // Copy results back to host + auto normalized_laplacian_data = std::vector(normalized_laplacian_csr_structure.get_nnz()); + auto normalized_laplacian_indices = + std::vector(normalized_laplacian_csr_structure.get_nnz()); + auto normalized_laplacian_indptr = + std::vector(normalized_laplacian_csr_structure.get_n_rows() + 1); + auto diagonal_data = std::vector(adjacency_structure_csr.get_n_rows()); + + raft::copy(&(normalized_laplacian_data[0]), + normalized_laplacian_csr.get_elements().data(), + normalized_laplacian_csr_structure.get_nnz(), + raft::resource::get_cuda_stream(res)); + raft::copy(&(normalized_laplacian_indices[0]), + normalized_laplacian_csr_structure.get_indices().data(), + normalized_laplacian_csr_structure.get_nnz(), + raft::resource::get_cuda_stream(res)); + raft::copy(&(normalized_laplacian_indptr[0]), + normalized_laplacian_csr_structure.get_indptr().data(), + normalized_laplacian_csr_structure.get_n_rows() + 1, + raft::resource::get_cuda_stream(res)); + raft::copy(&(diagonal_data[0]), + diagonal_out.data_handle(), + diagonal_out.size(), + raft::resource::get_cuda_stream(res)); + raft::resource::sync_stream(res); + + // For the given adjacency matrix, the diagonal degree matrix D has values [3, 2, 1, 2] + // The square root of these values is [√3, √2, 1, √2] + // The normalized Laplacian should have values close to: + // [1, -1/√(3*2), -1/√(3*1), -1/√(3*2), + // -1/√(2*3), 1, -1/√(2*2), + // -1/√(1*3), 1, + // -1/√(2*3), -1/√(2*2), 1] + + // Expected diagonal values (sqrt of diagonal degree matrix) + auto expected_diagonal = + std::vector{std::sqrt(3.0f), std::sqrt(2.0f), 1.0f, std::sqrt(2.0f)}; + + // Expected normalized Laplacian values + auto expected_data = std::vector{1.0f, + -1.0f / std::sqrt(3.0f * 2.0f), + -1.0f / std::sqrt(3.0f * 1.0f), + -1.0f / std::sqrt(3.0f * 2.0f), + -1.0f / std::sqrt(2.0f * 3.0f), + 1.0f, + -1.0f / std::sqrt(2.0f * 2.0f), + -1.0f / std::sqrt(1.0f * 3.0f), + 1.0f, + -1.0f / std::sqrt(2.0f * 3.0f), + -1.0f / std::sqrt(2.0f * 2.0f), + 1.0f}; + + // Same indices and indptr as non-normalized Laplacian + auto expected_indices = std::vector{0, 1, 2, 3, 0, 1, 3, 0, 2, 0, 1, 3}; + auto expected_indptr = std::vector{0, 4, 7, 9, 12}; + + // Compare results with expected values with a small tolerance for floating point differences + const float tol = 1e-6f; + ASSERT_EQ(expected_diagonal.size(), diagonal_data.size()); + for (size_t i = 0; i < expected_diagonal.size(); ++i) { + EXPECT_NEAR(expected_diagonal[i], diagonal_data[i], tol) + << "COO: Failed at diagonal index " << i; + } + + ASSERT_EQ(expected_data.size(), normalized_laplacian_data.size()); + for (size_t i = 0; i < expected_data.size(); ++i) { + EXPECT_NEAR(expected_data[i], normalized_laplacian_data[i], tol) + << "COO: Failed at data index " << i; + } + + EXPECT_EQ(expected_indices, normalized_laplacian_indices); + EXPECT_EQ(expected_indptr, normalized_laplacian_indptr); +} + } // namespace raft::sparse::linalg