From 68d412a119f918dc003ecd0c430b6bde7499bc96 Mon Sep 17 00:00:00 2001 From: William Hicks Date: Wed, 19 Feb 2025 00:13:33 -0500 Subject: [PATCH] Fix Laplacian calculation in spectral partitioning (#2568) These changes enable the use of the new Lanczos solver for spectral partitioning. This solver gives the correct eigenvalues in some cases where the old solver did not, but a previous attempt to introduce this led to a downstream breakage because the Laplacian was not being calculated correctly in the spectral partitioning logic. In this PR, we introduce a new `compute_graph_laplacian` function to correctly generate a new CSR matrix representing the graph Laplacian of the input CSR-formatted adjacency matrix. This resolves #2419. This function is used in the spectral partition function to obtain the correct partitioning using the new Lanczos solver. This eliminates all internal uses of `computeSmallestEigenvectors`, the old Lanczos implementation, which should allow it to be removed after a deprecation cycle for downstream consumers. This will in turn allow us to resolve #313. Note that the originally-reported breakage from the previous attempt to switch to the new Lanczos solver was used to create the test confirming correct functionality of the current PR. This PR is marked as a bugfix because it unblocks usage of the corrected Lanczos implementation. It is based on and requires PR #2541, and it will be marked as draft until that PR is merged. Authors: - William Hicks (https://github.com/wphicks) Approvers: - Victor Lafargue (https://github.com/viclafargue) - Divye Gala (https://github.com/divyegala) - Corey J. Nolet (https://github.com/cjnolet) URL: https://github.com/rapidsai/raft/pull/2568 --- cpp/include/raft/core/sparse_types.hpp | 2 +- .../raft/sparse/linalg/detail/laplacian.cuh | 118 ++++++++++++++++++ cpp/include/raft/sparse/linalg/laplacian.cuh | 39 ++++++ .../raft/spectral/detail/matrix_wrappers.hpp | 16 +++ .../raft/spectral/detail/partition.hpp | 8 +- cpp/include/raft/spectral/eigen_solvers.cuh | 53 +++++--- cpp/tests/CMakeLists.txt | 1 + cpp/tests/linalg/eigen_solvers.cu | 84 +++++++++++++ cpp/tests/sparse/laplacian.cu | 89 +++++++++++++ cpp/tests/sparse/spectral_matrix.cu | 2 + 10 files changed, 394 insertions(+), 18 deletions(-) create mode 100644 cpp/include/raft/sparse/linalg/detail/laplacian.cuh create mode 100644 cpp/include/raft/sparse/linalg/laplacian.cuh create mode 100644 cpp/tests/sparse/laplacian.cu diff --git a/cpp/include/raft/core/sparse_types.hpp b/cpp/include/raft/core/sparse_types.hpp index 6e5092f50f..8f3357c553 100644 --- a/cpp/include/raft/core/sparse_types.hpp +++ b/cpp/include/raft/core/sparse_types.hpp @@ -168,7 +168,7 @@ class sparse_matrix { row_type n_rows, col_type n_cols, nnz_type nnz = 0) noexcept(std::is_nothrow_default_constructible_v) - : structure_{handle, n_rows, n_cols, nnz}, cp_{}, c_elements_{cp_.create(handle, 0)} {}; + : structure_{handle, n_rows, n_cols, nnz}, cp_{}, c_elements_{cp_.create(handle, nnz)} {}; // Constructor that owns the data but not the structure // This constructor is only callable with a `structure_type == *_structure_view` diff --git a/cpp/include/raft/sparse/linalg/detail/laplacian.cuh b/cpp/include/raft/sparse/linalg/detail/laplacian.cuh new file mode 100644 index 0000000000..95318c74ed --- /dev/null +++ b/cpp/include/raft/sparse/linalg/detail/laplacian.cuh @@ -0,0 +1,118 @@ +/* + * Copyright (c) 2025, NVIDIA CORPORATION. + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +#pragma once +#include +#include +#include + +#include + +namespace raft { +namespace sparse { +namespace linalg { +namespace detail { + +/* Compute the graph Laplacian of an adjacency matrix + * + * This kernel implements the necessary logic for computing a graph + * Laplacian for an adjacency matrix in CSR format. A custom kernel is + * required because cusparse does not conveniently implement matrix subtraction with 64-bit + * indices. The custom kernel also allows the computation to be completed + * with no extra allocations or compute. + */ +template +RAFT_KERNEL compute_graph_laplacian_kernel(ElementType* output_values, + IndicesType* output_indices, + IndptrType* output_indptr, + IndptrType dim, + ElementType const* adj_values, + IndicesType const* adj_indices, + IndptrType const* adj_indptr) +{ + /* The graph Laplacian L of an adjacency matrix A is given by: + * L = D - A + * where D is the degree matrix of A. The degree matrix is itself defined + * as the sum of each row of A and represents the degree of the node + * indicated by the index of the row. */ + + for (auto row = threadIdx.x + blockIdx.x * blockDim.x; row < dim; row += blockDim.x * gridDim.x) { + auto row_begin = adj_indptr[row]; + auto row_end = adj_indptr[row + 1]; + // All output indexes will need to be offset by the row, since every row will + // gain exactly one new non-zero element. degree_output_index is the index + // where we will store the degree of each row + auto degree_output_index = row_begin + row; + auto degree_value = ElementType{}; + // value_index indicates the index of the current value in the original + // adjacency matrix + for (auto value_index = row_begin; value_index < row_end; ++value_index) { + auto col_index = adj_indices[value_index]; + auto is_lower_diagonal = col_index < row; + auto output_index = value_index + row + !is_lower_diagonal; + auto input_value = adj_values[value_index]; + degree_value += input_value; + output_values[output_index] = ElementType{-1} * input_value; + output_indices[output_index] = col_index; + // Increment the index where we will store the degree for every non-zero + // element before we reach the diagonal + degree_output_index += is_lower_diagonal; + } + output_values[degree_output_index] = degree_value; + output_indices[degree_output_index] = row; + output_indptr[row] = row_begin + row; + output_indptr[row + 1] = row_end + row + 1; + } +} + +template +auto compute_graph_laplacian( + raft::resources const& res, + device_csr_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 result = make_device_csr_matrix, + std::remove_const_t, + std::remove_const_t, + std::remove_const_t>( + res, + dim, + dim, + /* The nnz for the result will be the dimension of the (square) input matrix plus the number of + * non-zero elements in the original matrix, since we introduce non-zero elements along the + * diagonal to represent the degree of each node. */ + input_structure.get_nnz() + dim); + auto result_structure = result.structure_view(); + auto static constexpr const threads_per_block = 256; + auto blocks = std::min(int((dim + threads_per_block - 1) / threads_per_block), 65535); + auto stream = resource::get_cuda_stream(res); + detail::compute_graph_laplacian_kernel<<>>( + result.get_elements().data(), + result_structure.get_indices().data(), + result_structure.get_indptr().data(), + dim, + input.get_elements().data(), + input_structure.get_indices().data(), + input_structure.get_indptr().data()); + return result; +} + +} // namespace detail +} // namespace linalg +} // namespace sparse +} // namespace raft diff --git a/cpp/include/raft/sparse/linalg/laplacian.cuh b/cpp/include/raft/sparse/linalg/laplacian.cuh new file mode 100644 index 0000000000..69633f717c --- /dev/null +++ b/cpp/include/raft/sparse/linalg/laplacian.cuh @@ -0,0 +1,39 @@ +/* + * Copyright (c) 2025, NVIDIA CORPORATION. + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +#pragma once +#include +#include +#include + +namespace raft { +namespace sparse { +namespace linalg { + +/** Given a CSR adjacency matrix, return the graph Laplacian + * + * Note that for non-symmetric matrices, the out-degree Laplacian is returned. + */ +template +auto compute_graph_laplacian( + raft::resources const& res, + device_csr_matrix_view input) +{ + return detail::compute_graph_laplacian(res, input); +} + +} // namespace linalg +} // namespace sparse +} // namespace raft diff --git a/cpp/include/raft/spectral/detail/matrix_wrappers.hpp b/cpp/include/raft/spectral/detail/matrix_wrappers.hpp index af4a838bd5..417db8648d 100644 --- a/cpp/include/raft/spectral/detail/matrix_wrappers.hpp +++ b/cpp/include/raft/spectral/detail/matrix_wrappers.hpp @@ -15,6 +15,8 @@ */ #pragma once +#include +#include #include #include #include @@ -33,6 +35,7 @@ #include #include +#include // ========================================================= // Useful macros @@ -181,6 +184,19 @@ struct sparse_matrix_t { { } + auto to_csr_matrix_view() const + { + // The usage of sparse_matrix_t prior to introduction of this method + // assumed that all data was strictly on device. We will make the same + // assumption for construction of the csr_matrix_view + return device_csr_matrix_view{ + device_span{values_, std::uint64_t(nnz_)}, + device_compressed_structure_view{ + device_span{row_offsets_, std::uint64_t(nrows_ + 1)}, + device_span{col_indices_, std::uint64_t(nnz_)}, + ncols_}}; + } + virtual ~sparse_matrix_t(void) = default; // virtual because used as base for following matrix types diff --git a/cpp/include/raft/spectral/detail/partition.hpp b/cpp/include/raft/spectral/detail/partition.hpp index 26e4a73f9d..e99996275b 100644 --- a/cpp/include/raft/spectral/detail/partition.hpp +++ b/cpp/include/raft/spectral/detail/partition.hpp @@ -18,6 +18,7 @@ #include #include #include +#include #include #include #include @@ -97,14 +98,15 @@ std::tuple partition( // Compute eigenvectors of Laplacian // Initialize Laplacian - /// sparse_matrix_t A{handle, graph}; - spectral::matrix::laplacian_matrix_t L{handle, csr_m}; + auto laplacian = + raft::sparse::linalg::compute_graph_laplacian(handle, csr_m.to_csr_matrix_view()); auto eigen_config = eigen_solver.get_config(); auto nEigVecs = eigen_config.n_eigVecs; // Compute smallest eigenvalues and eigenvectors - std::get<0>(stats) = eigen_solver.solve_smallest_eigenvectors(handle, L, eigVals, eigVecs); + std::get<0>(stats) = + eigen_solver.solve_smallest_eigenvectors(handle, laplacian.view(), eigVals, eigVecs); // Whiten eigenvector matrix transform_eigen_matrix(handle, n, nEigVecs, eigVecs); diff --git a/cpp/include/raft/spectral/eigen_solvers.cuh b/cpp/include/raft/spectral/eigen_solvers.cuh index 03448e2b5e..0d511a2333 100644 --- a/cpp/include/raft/spectral/eigen_solvers.cuh +++ b/cpp/include/raft/spectral/eigen_solvers.cuh @@ -18,6 +18,7 @@ #pragma once +#include #include #include @@ -49,28 +50,52 @@ struct lanczos_solver_t { { } - index_type_t solve_smallest_eigenvectors( + auto solve_smallest_eigenvectors( raft::resources const& handle, matrix::sparse_matrix_t const& A, value_type_t* __restrict__ eigVals, value_type_t* __restrict__ eigVecs) const + { + auto csr_structure = + raft::make_device_compressed_structure_view( + const_cast(A.row_offsets_), + const_cast(A.col_indices_), + A.nrows_, + A.ncols_, + A.nnz_); + + auto csr_matrix = + raft::make_device_csr_matrix_view( + const_cast(A.values_), csr_structure); + return solve_smallest_eigenvectors(handle, csr_matrix, eigVals, eigVecs); + } + + auto solve_smallest_eigenvectors( + raft::resources const& handle, + device_csr_matrix_view input, + value_type_t* __restrict__ eigVals, + value_type_t* __restrict__ eigVecs) const { RAFT_EXPECTS(eigVals != nullptr, "Null eigVals buffer."); RAFT_EXPECTS(eigVecs != nullptr, "Null eigVecs buffer."); - index_type_t iters{}; - sparse::solver::computeSmallestEigenvectors(handle, - A, - config_.n_eigVecs, - config_.maxIter, - config_.restartIter, - config_.tol, - config_.reorthogonalize, - iters, - eigVals, - eigVecs, - config_.seed); - return iters; + auto lanczos_config = raft::sparse::solver::lanczos_solver_config{ + config_.n_eigVecs, config_.maxIter, config_.restartIter, config_.tol, config_.seed}; + auto v0_opt = std::optional>{ + std::nullopt}; + auto input_structure = input.structure_view(); + + auto solver_iterations = sparse::solver::lanczos_compute_smallest_eigenvectors( + handle, + lanczos_config, + input, + v0_opt, + raft::make_device_vector_view(eigVals, + config_.n_eigVecs), + raft::make_device_matrix_view( + eigVecs, input_structure.get_n_rows(), config_.n_eigVecs)); + + return index_type_t{solver_iterations}; } index_type_t solve_largest_eigenvectors( diff --git a/cpp/tests/CMakeLists.txt b/cpp/tests/CMakeLists.txt index 34bea67fbe..12da3bec13 100644 --- a/cpp/tests/CMakeLists.txt +++ b/cpp/tests/CMakeLists.txt @@ -259,6 +259,7 @@ if(BUILD_TESTS) sparse/csr_transpose.cu sparse/degree.cu sparse/filter.cu + sparse/laplacian.cu sparse/masked_matmul.cu sparse/norm.cu sparse/normalize.cu diff --git a/cpp/tests/linalg/eigen_solvers.cu b/cpp/tests/linalg/eigen_solvers.cu index 250deb6f33..f43081690c 100644 --- a/cpp/tests/linalg/eigen_solvers.cu +++ b/cpp/tests/linalg/eigen_solvers.cu @@ -14,12 +14,16 @@ * limitations under the License. */ +#include #include #include #include #include +#include #include +#include + #include #include @@ -117,5 +121,85 @@ TEST(Raft, SpectralSolvers) EXPECT_ANY_THROW(spectral::analyzePartition(h, sm, k, clusters, edgeCut, cost)); } +TEST(Raft, SpectralPartition) +{ + auto offsets = thrust::device_vector(std::vector{ + 0, 16, 25, 35, 41, 44, 48, 52, 56, 61, 63, 66, 67, 69, 74, 76, 78, 80, + 82, 84, 87, 89, 91, 93, 98, 101, 104, 106, 110, 113, 117, 121, 127, 139, 156}); + auto indices = thrust::device_vector(std::vector{ + 1, 2, 3, 4, 5, 6, 7, 8, 10, 11, 12, 13, 17, 19, 21, 31, 0, 2, 3, 7, 13, 17, 19, + 21, 30, 0, 1, 3, 7, 8, 9, 13, 27, 28, 32, 0, 1, 2, 7, 12, 13, 0, 6, 10, 0, 6, + 10, 16, 0, 4, 5, 16, 0, 1, 2, 3, 0, 2, 30, 32, 33, 2, 33, 0, 4, 5, 0, 0, 3, + 0, 1, 2, 3, 33, 32, 33, 32, 33, 5, 6, 0, 1, 32, 33, 0, 1, 33, 32, 33, 0, 1, 32, + 33, 25, 27, 29, 32, 33, 25, 27, 31, 23, 24, 31, 29, 33, 2, 23, 24, 33, 2, 31, 33, 23, 26, + 32, 33, 1, 8, 32, 33, 0, 24, 25, 28, 32, 33, 2, 8, 14, 15, 18, 20, 22, 23, 29, 30, 31, + 33, 8, 9, 13, 14, 15, 18, 19, 20, 22, 23, 26, 27, 28, 29, 30, 31, 32}); + auto values = thrust::device_vector(std::vector{ + 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, + 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, + 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, + 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, + 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, + 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, + 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, + 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, + 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0}); + + auto num_verts = int(offsets.size() - 1); + auto num_edges = indices.size(); + + auto result_v = thrust::device_vector(std::vector(num_verts, -1)); + + auto constexpr const n_clusters = 8; + auto constexpr const n_eig_vects = std::uint64_t{8}; + + auto constexpr const evs_tolerance = .00001f; + auto constexpr const kmean_tolerance = .00001f; + auto constexpr const evs_max_iter = 100; + auto constexpr const kmean_max_iter = 100; + + auto eig_vals = thrust::device_vector(n_eig_vects); + auto eig_vects = thrust::device_vector(n_eig_vects * num_verts); + + auto handle = raft::handle_t{}; + + auto restartIter_lanczos = int{15 + n_eig_vects}; + + auto seed_eig_solver = std::uint64_t{1234567}; + auto seed_cluster_solver = std::uint64_t{12345678}; + + auto const adj_matrix = raft::spectral::matrix::sparse_matrix_t{handle, + offsets.data().get(), + indices.data().get(), + values.data().get(), + num_verts, + num_verts, + num_edges}; + + auto eig_cfg = raft::spectral::eigen_solver_config_t{ + n_eig_vects, evs_max_iter, restartIter_lanczos, evs_tolerance, false, seed_eig_solver}; + auto eigen_solver = raft::spectral::lanczos_solver_t{eig_cfg}; + + auto clust_cfg = raft::spectral::cluster_solver_config_t{ + n_clusters, kmean_max_iter, kmean_tolerance, seed_cluster_solver}; + auto cluster_solver = raft::spectral::kmeans_solver_t{clust_cfg}; + + partition(handle, + adj_matrix, + eigen_solver, + cluster_solver, + result_v.data().get(), + eig_vals.data().get(), + eig_vects.data().get()); + + auto edge_cut = float{}; + auto cost = float{}; + + raft::spectral::analyzePartition( + handle, adj_matrix, n_clusters, result_v.data().get(), edge_cut, cost); + + ASSERT_LT(edge_cut, 55.0); +} + } // namespace spectral } // namespace raft diff --git a/cpp/tests/sparse/laplacian.cu b/cpp/tests/sparse/laplacian.cu new file mode 100644 index 0000000000..ea1d158f3e --- /dev/null +++ b/cpp/tests/sparse/laplacian.cu @@ -0,0 +1,89 @@ +/* + * Copyright (c) 2020-2024, NVIDIA CORPORATION. + * + * Licensed under the Apache License, Version 2.0 (the "License"); + * you may not use this file except in compliance with the License. + * You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ + +#include +#include +#include +#include +#include + +#include + +#include + +namespace raft { +namespace sparse { +namespace linalg { + +TEST(Raft, ComputeGraphLaplacian) +{ + // The following adjacency matrix will be used to allow for manual + // verification of results: + // [[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{}; + auto adjacency_matrix = + make_device_csr_matrix(res, int(indptr.size() - 1), int(indptr.size() - 1), data.size()); + auto adjacency_structure = adjacency_matrix.structure_view(); + raft::copy(adjacency_matrix.get_elements().data(), + &(data[0]), + data.size(), + raft::resource::get_cuda_stream(res)); + raft::copy(adjacency_structure.get_indices().data(), + &(indices[0]), + indices.size(), + raft::resource::get_cuda_stream(res)); + raft::copy(adjacency_structure.get_indptr().data(), + &(indptr[0]), + indptr.size(), + raft::resource::get_cuda_stream(res)); + auto laplacian = compute_graph_laplacian(res, adjacency_matrix.view()); + auto laplacian_structure = laplacian.structure_view(); + auto laplacian_data = std::vector(laplacian_structure.get_nnz()); + auto laplacian_indices = std::vector(laplacian_structure.get_nnz()); + auto laplacian_indptr = std::vector(laplacian_structure.get_n_rows() + 1); + raft::copy(&(laplacian_data[0]), + laplacian.get_elements().data(), + laplacian_structure.get_nnz(), + raft::resource::get_cuda_stream(res)); + raft::copy(&(laplacian_indices[0]), + laplacian_structure.get_indices().data(), + laplacian_structure.get_nnz(), + raft::resource::get_cuda_stream(res)); + raft::copy(&(laplacian_indptr[0]), + laplacian_structure.get_indptr().data(), + laplacian_structure.get_n_rows() + 1, + raft::resource::get_cuda_stream(res)); + auto expected_data = std::vector{3, -1, -1, -1, -1, 2, -1, -1, 1, -1, -1, 2}; + 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}; + raft::resource::sync_stream(res); + + EXPECT_EQ(expected_data, laplacian_data); + EXPECT_EQ(expected_indices, laplacian_indices); + EXPECT_EQ(expected_indptr, laplacian_indptr); +} + +} // namespace linalg +} // namespace sparse +} // namespace raft diff --git a/cpp/tests/sparse/spectral_matrix.cu b/cpp/tests/sparse/spectral_matrix.cu index 6d52dfb1bb..962fd25d04 100644 --- a/cpp/tests/sparse/spectral_matrix.cu +++ b/cpp/tests/sparse/spectral_matrix.cu @@ -17,7 +17,9 @@ #include #include #include +#include #include +#include #include