diff --git a/cpp/CMakeLists.txt b/cpp/CMakeLists.txt index d4e9f2c6b4..c138eae53b 100644 --- a/cpp/CMakeLists.txt +++ b/cpp/CMakeLists.txt @@ -664,6 +664,7 @@ if(NOT BUILD_CPU_ONLY) src/preprocessing/quantize/binary.cu src/preprocessing/quantize/pq.cu src/preprocessing/spectral/spectral_embedding.cu + src/preprocessing/pca/pca.cu src/selection/select_k_float_int64_t.cu src/selection/select_k_float_int32_t.cu src/selection/select_k_float_uint32_t.cu diff --git a/cpp/include/cuvs/preprocessing/pca.hpp b/cpp/include/cuvs/preprocessing/pca.hpp new file mode 100644 index 0000000000..f30b245d34 --- /dev/null +++ b/cpp/include/cuvs/preprocessing/pca.hpp @@ -0,0 +1,185 @@ +/* + * SPDX-FileCopyrightText: Copyright (c) 2026, NVIDIA CORPORATION. + * SPDX-License-Identifier: Apache-2.0 + */ + +#pragma once + +#include +#include +#include + +namespace cuvs::preprocessing::pca { + +using solver = raft::linalg::solver; + +/** + * @brief Parameters for PCA decomposition. Ref: + * http://scikit-learn.org/stable/modules/generated/sklearn.decomposition.PCA.html + */ +struct params { + /** @brief Number of components to keep. */ + int n_components = 1; + + /** + * @brief If false, data passed to fit are overwritten and running fit(X).transform(X) will + * not yield the expected results, use fit_transform(X) instead. + */ + bool copy = true; + + /** + * @brief When true (false by default) the components vectors are multiplied by the square + * root of n_samples and then divided by the singular values to ensure uncorrelated outputs with + * unit component-wise variances. + */ + bool whiten = false; + + /** @brief The solver algorithm to use. */ + solver algorithm = solver::COV_EIG_DQ; + + /** + * @brief Tolerance for singular values computed by svd_solver == 'arpack' or + * the Jacobi solver. + */ + float tol = 0.0f; + + /** + * @brief Number of iterations for the power method computed by the Jacobi solver. + */ + int n_iterations = 15; +}; + +/** + * @defgroup pca PCA (Principal Component Analysis) + * @{ + */ + +/** + * @brief Perform PCA fit operation. + * + * Computes the principal components, explained variances, singular values, and column means + * from the input data. + * + * @code{.cpp} + * #include + * #include + * + * raft::resources handle; + * + * cuvs::preprocessing::pca::params params; + * params.n_components = 2; + * + * auto input = raft::make_device_matrix(handle, n_rows, n_cols); + * // ... fill input ... + * + * auto components = raft::make_device_matrix( + * handle, params.n_components, n_cols); + * auto explained_var = raft::make_device_vector(handle, params.n_components); + * auto explained_var_ratio = raft::make_device_vector(handle, params.n_components); + * auto singular_vals = raft::make_device_vector(handle, params.n_components); + * auto mu = raft::make_device_vector(handle, n_cols); + * auto noise_vars = raft::make_device_scalar(handle); + * + * cuvs::preprocessing::pca::fit(handle, params, + * input.view(), components.view(), explained_var.view(), + * explained_var_ratio.view(), singular_vals.view(), mu.view(), noise_vars.view()); + * @endcode + * + * @param[in] handle raft resource handle + * @param[in] config PCA parameters + * @param[inout] input input data [n_rows x n_cols] (col-major). Modified temporarily. + * @param[out] components principal components [n_components x n_cols] (col-major) + * @param[out] explained_var explained variances [n_components] + * @param[out] explained_var_ratio explained variance ratios [n_components] + * @param[out] singular_vals singular values [n_components] + * @param[out] mu column means [n_cols] + * @param[out] noise_vars noise variance (scalar) + * @param[in] flip_signs_based_on_U whether to determine signs by U (true) or V.T (false) + */ +void fit(raft::resources const& handle, + const params& config, + raft::device_matrix_view input, + raft::device_matrix_view components, + raft::device_vector_view explained_var, + raft::device_vector_view explained_var_ratio, + raft::device_vector_view singular_vals, + raft::device_vector_view mu, + raft::device_scalar_view noise_vars, + bool flip_signs_based_on_U = false); + +/** + * @brief Perform PCA fit and transform operations. + * + * Computes the principal components and transforms the input data into the eigenspace + * in a single operation. + * + * @param[in] handle raft resource handle + * @param[in] config PCA parameters + * @param[inout] input input data [n_rows x n_cols] (col-major). Modified temporarily. + * @param[out] trans_input transformed data [n_rows x n_components] (col-major) + * @param[out] components principal components [n_components x n_cols] (col-major) + * @param[out] explained_var explained variances [n_components] + * @param[out] explained_var_ratio explained variance ratios [n_components] + * @param[out] singular_vals singular values [n_components] + * @param[out] mu column means [n_cols] + * @param[out] noise_vars noise variance (scalar) + * @param[in] flip_signs_based_on_U whether to determine signs by U (true) or V.T (false) + */ +void fit_transform(raft::resources const& handle, + const params& config, + raft::device_matrix_view input, + raft::device_matrix_view trans_input, + raft::device_matrix_view components, + raft::device_vector_view explained_var, + raft::device_vector_view explained_var_ratio, + raft::device_vector_view singular_vals, + raft::device_vector_view mu, + raft::device_scalar_view noise_vars, + bool flip_signs_based_on_U = false); + +/** + * @brief Perform PCA transform operation. + * + * Transforms the input data into the eigenspace using previously computed principal components. + * + * @param[in] handle raft resource handle + * @param[in] config PCA parameters + * @param[inout] input data to transform [n_rows x n_cols] (col-major). Modified temporarily + * (mean-centered then restored). + * @param[in] components principal components [n_components x n_cols] (col-major) + * @param[in] singular_vals singular values [n_components] + * @param[in] mu column means [n_cols] + * @param[out] trans_input transformed data [n_rows x n_components] (col-major) + */ +void transform(raft::resources const& handle, + const params& config, + raft::device_matrix_view input, + raft::device_matrix_view components, + raft::device_vector_view singular_vals, + raft::device_vector_view mu, + raft::device_matrix_view trans_input); + +/** + * @brief Perform PCA inverse transform operation. + * + * Transforms data from the eigenspace back to the original space. + * + * @param[in] handle raft resource handle + * @param[in] config PCA parameters + * @param[in] trans_input transformed data [n_rows x n_components] (col-major) + * @param[in] components principal components [n_components x n_cols] (col-major) + * @param[in] singular_vals singular values [n_components] + * @param[in] mu column means [n_cols] + * @param[out] output reconstructed data [n_rows x n_cols] (col-major) + */ +void inverse_transform(raft::resources const& handle, + const params& config, + raft::device_matrix_view trans_input, + raft::device_matrix_view components, + raft::device_vector_view singular_vals, + raft::device_vector_view mu, + raft::device_matrix_view output); + +/** @} */ // end group pca + +} // namespace cuvs::preprocessing::pca diff --git a/cpp/src/preprocessing/pca/detail/pca.cuh b/cpp/src/preprocessing/pca/detail/pca.cuh new file mode 100644 index 0000000000..0aa1f555e2 --- /dev/null +++ b/cpp/src/preprocessing/pca/detail/pca.cuh @@ -0,0 +1,113 @@ +/* + * SPDX-FileCopyrightText: Copyright (c) 2026, NVIDIA CORPORATION. + * SPDX-License-Identifier: Apache-2.0 + */ + +#pragma once + +#include + +#include +#include +#include + +namespace cuvs::preprocessing::pca::detail { + +/** + * @brief Convert cuvs::preprocessing::pca::params to raft::linalg::paramsPCA. + */ +inline auto to_raft_params(const params& config, std::size_t n_rows, std::size_t n_cols) + -> raft::linalg::paramsPCA +{ + raft::linalg::paramsPCA prms; + prms.n_rows = n_rows; + prms.n_cols = n_cols; + prms.n_components = config.n_components; + prms.algorithm = config.algorithm; + prms.tol = config.tol; + prms.n_iterations = config.n_iterations; + prms.copy = config.copy; + prms.whiten = config.whiten; + return prms; +} + +template +void fit(raft::resources const& handle, + const params& config, + raft::device_matrix_view input, + raft::device_matrix_view components, + raft::device_vector_view explained_var, + raft::device_vector_view explained_var_ratio, + raft::device_vector_view singular_vals, + raft::device_vector_view mu, + raft::device_scalar_view noise_vars, + bool flip_signs_based_on_U) +{ + auto raft_prms = to_raft_params(config, input.extent(0), input.extent(1)); + raft::linalg::pca_fit(handle, + raft_prms, + input, + components, + explained_var, + explained_var_ratio, + singular_vals, + mu, + noise_vars, + flip_signs_based_on_U); +} + +template +void fit_transform(raft::resources const& handle, + const params& config, + raft::device_matrix_view input, + raft::device_matrix_view trans_input, + raft::device_matrix_view components, + raft::device_vector_view explained_var, + raft::device_vector_view explained_var_ratio, + raft::device_vector_view singular_vals, + raft::device_vector_view mu, + raft::device_scalar_view noise_vars, + bool flip_signs_based_on_U) +{ + auto raft_prms = to_raft_params(config, input.extent(0), input.extent(1)); + raft::linalg::pca_fit_transform(handle, + raft_prms, + input, + trans_input, + components, + explained_var, + explained_var_ratio, + singular_vals, + mu, + noise_vars, + flip_signs_based_on_U); +} + +template +void transform(raft::resources const& handle, + const params& config, + raft::device_matrix_view input, + raft::device_matrix_view components, + raft::device_vector_view singular_vals, + raft::device_vector_view mu, + raft::device_matrix_view trans_input) +{ + auto raft_prms = to_raft_params(config, input.extent(0), input.extent(1)); + raft::linalg::pca_transform(handle, raft_prms, input, components, singular_vals, mu, trans_input); +} + +template +void inverse_transform(raft::resources const& handle, + const params& config, + raft::device_matrix_view trans_input, + raft::device_matrix_view components, + raft::device_vector_view singular_vals, + raft::device_vector_view mu, + raft::device_matrix_view output) +{ + auto raft_prms = to_raft_params(config, output.extent(0), output.extent(1)); + raft::linalg::pca_inverse_transform( + handle, raft_prms, trans_input, components, singular_vals, mu, output); +} + +} // namespace cuvs::preprocessing::pca::detail diff --git a/cpp/src/preprocessing/pca/pca.cu b/cpp/src/preprocessing/pca/pca.cu new file mode 100644 index 0000000000..ac944fddd7 --- /dev/null +++ b/cpp/src/preprocessing/pca/pca.cu @@ -0,0 +1,98 @@ +/* + * SPDX-FileCopyrightText: Copyright (c) 2026, NVIDIA CORPORATION. + * SPDX-License-Identifier: Apache-2.0 + */ + +#include "./detail/pca.cuh" + +#include + +namespace cuvs::preprocessing::pca { + +#define CUVS_INST_PCA_FIT(DataT, IndexT) \ + void fit(raft::resources const& handle, \ + const params& config, \ + raft::device_matrix_view input, \ + raft::device_matrix_view components, \ + raft::device_vector_view explained_var, \ + raft::device_vector_view explained_var_ratio, \ + raft::device_vector_view singular_vals, \ + raft::device_vector_view mu, \ + raft::device_scalar_view noise_vars, \ + bool flip_signs_based_on_U) \ + { \ + detail::fit(handle, \ + config, \ + input, \ + components, \ + explained_var, \ + explained_var_ratio, \ + singular_vals, \ + mu, \ + noise_vars, \ + flip_signs_based_on_U); \ + } + +CUVS_INST_PCA_FIT(float, int64_t); +#undef CUVS_INST_PCA_FIT + +#define CUVS_INST_PCA_FIT_TRANSFORM(DataT, IndexT) \ + void fit_transform(raft::resources const& handle, \ + const params& config, \ + raft::device_matrix_view input, \ + raft::device_matrix_view trans_input, \ + raft::device_matrix_view components, \ + raft::device_vector_view explained_var, \ + raft::device_vector_view explained_var_ratio, \ + raft::device_vector_view singular_vals, \ + raft::device_vector_view mu, \ + raft::device_scalar_view noise_vars, \ + bool flip_signs_based_on_U) \ + { \ + detail::fit_transform(handle, \ + config, \ + input, \ + trans_input, \ + components, \ + explained_var, \ + explained_var_ratio, \ + singular_vals, \ + mu, \ + noise_vars, \ + flip_signs_based_on_U); \ + } + +CUVS_INST_PCA_FIT_TRANSFORM(float, int64_t); +#undef CUVS_INST_PCA_FIT_TRANSFORM + +#define CUVS_INST_PCA_TRANSFORM(DataT, IndexT) \ + void transform(raft::resources const& handle, \ + const params& config, \ + raft::device_matrix_view input, \ + raft::device_matrix_view components, \ + raft::device_vector_view singular_vals, \ + raft::device_vector_view mu, \ + raft::device_matrix_view trans_input) \ + { \ + detail::transform(handle, config, input, components, singular_vals, mu, trans_input); \ + } + +CUVS_INST_PCA_TRANSFORM(float, int64_t); +#undef CUVS_INST_PCA_TRANSFORM + +#define CUVS_INST_PCA_INVERSE_TRANSFORM(DataT, IndexT) \ + void inverse_transform(raft::resources const& handle, \ + const params& config, \ + raft::device_matrix_view trans_input, \ + raft::device_matrix_view components, \ + raft::device_vector_view singular_vals, \ + raft::device_vector_view mu, \ + raft::device_matrix_view output) \ + { \ + detail::inverse_transform(handle, config, trans_input, components, singular_vals, mu, output); \ + } + +CUVS_INST_PCA_INVERSE_TRANSFORM(float, int64_t); +#undef CUVS_INST_PCA_INVERSE_TRANSFORM + +} // namespace cuvs::preprocessing::pca diff --git a/cpp/tests/CMakeLists.txt b/cpp/tests/CMakeLists.txt index bbddef87e5..62fbf960ab 100644 --- a/cpp/tests/CMakeLists.txt +++ b/cpp/tests/CMakeLists.txt @@ -339,8 +339,11 @@ ConfigureTest( ConfigureTest( NAME PREPROCESSING_TEST - PATH preprocessing/scalar_quantization.cu preprocessing/binary_quantization.cu - preprocessing/spectral_embedding.cu preprocessing/product_quantization.cu + PATH preprocessing/scalar_quantization.cu + preprocessing/binary_quantization.cu + preprocessing/spectral_embedding.cu + preprocessing/product_quantization.cu + preprocessing/pca.cu GPUS 1 PERCENT 100 ) diff --git a/cpp/tests/preprocessing/pca.cu b/cpp/tests/preprocessing/pca.cu new file mode 100644 index 0000000000..4430972911 --- /dev/null +++ b/cpp/tests/preprocessing/pca.cu @@ -0,0 +1,266 @@ +/* + * SPDX-FileCopyrightText: Copyright (c) 2018-2026, NVIDIA CORPORATION. + * SPDX-License-Identifier: Apache-2.0 + */ + +#include "../test_utils.cuh" + +#include + +#include +#include +#include +#include +#include + +#include +#include + +#include +#include +#include + +namespace cuvs::preprocessing::pca { + +template +struct PcaInputs { + T tolerance; + int n_row; + int n_col; + int n_row2; + int n_col2; + unsigned long long int seed; + int algo; +}; + +template +::std::ostream& operator<<(::std::ostream& os, const PcaInputs& dims) +{ + return os; +} + +/** + * @brief Run fit_transform followed by inverse_transform. + * + * Intermediate buffers are managed internally unless the caller provides + * pre-allocated pointers via the optional parameters, in which case the + * results are written there directly. + */ +template +void pca_roundtrip(raft::resources const& handle, + T* input, + int n_rows, + int n_cols, + T* output, + int n_components, + int algo, + cudaStream_t stream, + T* components_out = nullptr, + T* explained_var_out = nullptr, + T* trans_out = nullptr) +{ + params prms; + prms.n_components = n_components; + if (algo == 0) + prms.algorithm = solver::COV_EIG_DQ; + else + prms.algorithm = solver::COV_EIG_JACOBI; + + rmm::device_uvector comp_buf(components_out ? 0 : n_components * n_cols, stream); + rmm::device_uvector ev_buf(explained_var_out ? 0 : n_components, stream); + rmm::device_uvector trans_buf(trans_out ? 0 : n_rows * n_components, stream); + + T* comp_ptr = components_out ? components_out : comp_buf.data(); + T* ev_ptr = explained_var_out ? explained_var_out : ev_buf.data(); + T* trans_ptr = trans_out ? trans_out : trans_buf.data(); + + rmm::device_uvector evr(n_components, stream); + rmm::device_uvector sv(n_components, stream); + rmm::device_uvector mu(n_cols, stream); + rmm::device_uvector nv(1, stream); + + auto input_view = + raft::make_device_matrix_view(input, n_rows, n_cols); + auto trans_view = + raft::make_device_matrix_view(trans_ptr, n_rows, n_components); + auto comp_view = + raft::make_device_matrix_view(comp_ptr, n_components, n_cols); + auto ev_view = raft::make_device_vector_view(ev_ptr, n_components); + auto evr_view = raft::make_device_vector_view(evr.data(), n_components); + auto sv_view = raft::make_device_vector_view(sv.data(), n_components); + auto mu_view = raft::make_device_vector_view(mu.data(), n_cols); + auto nv_view = raft::make_device_scalar_view(nv.data()); + auto output_view = + raft::make_device_matrix_view(output, n_rows, n_cols); + + fit_transform( + handle, prms, input_view, trans_view, comp_view, ev_view, evr_view, sv_view, mu_view, nv_view); + inverse_transform(handle, prms, trans_view, comp_view, sv_view, mu_view, output_view); +} + +template +class PcaTest : public ::testing::TestWithParam> { + public: + PcaTest() + : params_(::testing::TestWithParam>::GetParam()), + stream(raft::resource::get_cuda_stream(handle)), + explained_vars(params_.n_col, stream), + explained_vars_ref(params_.n_col, stream), + components(params_.n_col * params_.n_col, stream), + components_ref(params_.n_col * params_.n_col, stream), + trans_data(params_.n_row * params_.n_col, stream), + trans_data_ref(params_.n_row * params_.n_col, stream), + data(params_.n_row * params_.n_col, stream), + data_back(params_.n_row * params_.n_col, stream), + data2(params_.n_row2 * params_.n_col2, stream), + data2_back(params_.n_row2 * params_.n_col2, stream) + { + } + + protected: + void SetUp() override + { + int len = params_.n_row * params_.n_col; + int len2 = params_.n_row2 * params_.n_col2; + + // --- basic test: all components, known reference data --- + { + std::vector data_h = {1.0, 2.0, 5.0, 4.0, 2.0, 1.0}; + data_h.resize(len); + raft::update_device(data.data(), data_h.data(), len, stream); + + std::vector trans_data_ref_h = {-2.3231, -0.3517, 2.6748, 0.3979, -0.6571, 0.2592}; + trans_data_ref_h.resize(len); + raft::update_device(trans_data_ref.data(), trans_data_ref_h.data(), len, stream); + + int len_comp = params_.n_col * params_.n_col; + + std::vector components_ref_h = {0.8163, 0.5776, -0.5776, 0.8163}; + components_ref_h.resize(len_comp); + std::vector explained_vars_ref_h = {6.338, 0.3287}; + explained_vars_ref_h.resize(params_.n_col); + + raft::update_device(components_ref.data(), components_ref_h.data(), len_comp, stream); + raft::update_device( + explained_vars_ref.data(), explained_vars_ref_h.data(), params_.n_col, stream); + + pca_roundtrip(handle, + data.data(), + params_.n_row, + params_.n_col, + data_back.data(), + params_.n_col, + params_.algo, + stream, + components.data(), + explained_vars.data(), + trans_data.data()); + } + + // --- advanced test: all components, random data --- + { + raft::random::Rng r(params_.seed, raft::random::GenPC); + r.uniform(data2.data(), len2, T(-1.0), T(1.0), stream); + + pca_roundtrip(handle, + data2.data(), + params_.n_row2, + params_.n_col2, + data2_back.data(), + params_.n_col2, + params_.algo, + stream); + } + + // --- dim reduction test: n_components < n_cols, random data --- + { + int n_components = std::max(1, params_.n_col2 / 4); + + rmm::device_uvector input(len2, stream); + rmm::device_uvector input_copy(len2, stream); + rmm::device_uvector recon(len2, stream); + + raft::random::Rng rng(params_.seed + 1, raft::random::GenPC); + rng.uniform(input.data(), len2, T(-1.0), T(1.0), stream); + raft::copy(input_copy.data(), input.data(), len2, stream); + + pca_roundtrip(handle, + input.data(), + params_.n_row2, + params_.n_col2, + recon.data(), + n_components, + params_.algo, + stream); + + std::vector orig_h(len2); + std::vector recon_h(len2); + raft::update_host(orig_h.data(), input_copy.data(), len2, stream); + raft::update_host(recon_h.data(), recon.data(), len2, stream); + RAFT_CUDA_TRY(cudaStreamSynchronize(stream)); + + max_recon_err = T(0); + for (int i = 0; i < len2; ++i) { + max_recon_err = std::max(max_recon_err, std::abs(orig_h[i] - recon_h[i])); + } + } + } + + void testPca() + { + auto s = raft::resource::get_cuda_stream(handle); + + ASSERT_TRUE(devArrMatch(explained_vars.data(), + explained_vars_ref.data(), + params_.n_col, + cuvs::CompareApprox(params_.tolerance), + s)); + + ASSERT_TRUE(devArrMatch(components.data(), + components_ref.data(), + (params_.n_col * params_.n_col), + cuvs::CompareApprox(params_.tolerance), + s)); + + ASSERT_TRUE(devArrMatch(trans_data.data(), + trans_data_ref.data(), + (params_.n_row * params_.n_col), + cuvs::CompareApprox(params_.tolerance), + s)); + + ASSERT_TRUE(devArrMatch(data.data(), + data_back.data(), + (params_.n_row * params_.n_col), + cuvs::CompareApprox(params_.tolerance), + s)); + + ASSERT_TRUE(devArrMatch(data2.data(), + data2_back.data(), + (params_.n_row2 * params_.n_col2), + cuvs::CompareApprox(params_.tolerance), + s)); + + EXPECT_GT(max_recon_err, T(1e-5)) << "Error should be non-zero when n_components < n_cols"; + EXPECT_LT(max_recon_err, T(2.0)) << "Reconstruction error should be bounded"; + } + + private: + raft::device_resources handle; + cudaStream_t stream; + + PcaInputs params_; + T max_recon_err = T(0); + + rmm::device_uvector explained_vars, explained_vars_ref, components, components_ref, trans_data, + trans_data_ref, data, data_back, data2, data2_back; +}; + +const std::vector> inputsf2 = {{0.01f, 3, 2, 1024, 128, 1234ULL, 0}, + {0.01f, 3, 2, 256, 32, 1234ULL, 1}}; + +typedef PcaTest PcaTestF; +TEST_P(PcaTestF, Result) { this->testPca(); } + +INSTANTIATE_TEST_CASE_P(PcaTests, PcaTestF, ::testing::ValuesIn(inputsf2)); + +} // end namespace cuvs::preprocessing::pca