diff --git a/cpp/include/cuml/decomposition/params.hpp b/cpp/include/cuml/decomposition/params.hpp index edbf005a05..aec8e7a67a 100644 --- a/cpp/include/cuml/decomposition/params.hpp +++ b/cpp/include/cuml/decomposition/params.hpp @@ -5,20 +5,14 @@ #pragma once +#include + #include namespace ML { -/** - * @param COV_EIG_DQ: covariance of input will be used along with eigen decomposition using divide - * and conquer method for symmetric matrices - * @param COV_EIG_JACOBI: covariance of input will be used along with eigen decomposition using - * jacobi method for symmetric matrices - */ -enum class solver : int { - COV_EIG_DQ, - COV_EIG_JACOBI, -}; +using solver = raft::linalg::solver; +using mg_solver = raft::linalg::solver; class params { public: @@ -53,10 +47,8 @@ class paramsTSVDTemplate : public paramsSolver { * root of n_samples and then divided by the singular values to ensure uncorrelated outputs with * unit component-wise variances. * @param algorithm: the solver to be used in PCA. - * @param tol: Tolerance for singular values computed by svd_solver == ‘arpack’ or svd_solver == - * ‘COV_EIG_JACOBI’ - * @param n_iterations: Number of iterations for the power method computed by jacobi method - * (svd_solver == 'COV_EIG_JACOBI'). + * @param tol: Tolerance for singular values computed by the Jacobi solver + * @param n_iterations: Number of iterations for the power method computed by the Jacobi solver * @param verbose: 0: no error message printing, 1: print error messages */ @@ -70,9 +62,29 @@ class paramsPCATemplate : public paramsTSVDTemplate { typedef paramsTSVDTemplate<> paramsTSVD; typedef paramsPCATemplate<> paramsPCA; -enum class mg_solver { COV_EIG_DQ, COV_EIG_JACOBI }; - typedef paramsPCATemplate paramsPCAMG; typedef paramsTSVDTemplate paramsTSVDMG; +template +inline raft::linalg::paramsTSVD to_raft_params(const paramsTSVDTemplate& ml_prms) +{ + raft::linalg::paramsTSVD prms; + prms.tol = ml_prms.tol; + prms.n_iterations = ml_prms.n_iterations; + prms.algorithm = ml_prms.algorithm; + return prms; +} + +template +inline raft::linalg::paramsPCA to_raft_params(const paramsPCATemplate& ml_prms) +{ + raft::linalg::paramsPCA prms; + prms.tol = ml_prms.tol; + prms.n_iterations = ml_prms.n_iterations; + prms.algorithm = ml_prms.algorithm; + prms.copy = ml_prms.copy; + prms.whiten = ml_prms.whiten; + return prms; +} + }; // end namespace ML diff --git a/cpp/include/cuml/decomposition/pca.hpp b/cpp/include/cuml/decomposition/pca.hpp index 875617387f..1c29f44d0a 100644 --- a/cpp/include/cuml/decomposition/pca.hpp +++ b/cpp/include/cuml/decomposition/pca.hpp @@ -1,5 +1,5 @@ /* - * SPDX-FileCopyrightText: Copyright (c) 2018-2025, NVIDIA CORPORATION. + * SPDX-FileCopyrightText: Copyright (c) 2018-2026, NVIDIA CORPORATION. * SPDX-License-Identifier: Apache-2.0 */ @@ -13,7 +13,7 @@ class handle_t; namespace ML { -void pcaFit(raft::handle_t& handle, +void pcaFit(const raft::handle_t& handle, float* input, float* components, float* explained_var, @@ -23,7 +23,7 @@ void pcaFit(raft::handle_t& handle, float* noise_vars, const paramsPCA& prms, bool flip_signs_based_on_U); -void pcaFit(raft::handle_t& handle, +void pcaFit(const raft::handle_t& handle, double* input, double* components, double* explained_var, @@ -33,7 +33,7 @@ void pcaFit(raft::handle_t& handle, double* noise_vars, const paramsPCA& prms, bool flip_signs_based_on_U); -void pcaFitTransform(raft::handle_t& handle, +void pcaFitTransform(const raft::handle_t& handle, float* input, float* trans_input, float* components, @@ -44,7 +44,7 @@ void pcaFitTransform(raft::handle_t& handle, float* noise_vars, const paramsPCA& prms, bool flip_signs_based_on_U); -void pcaFitTransform(raft::handle_t& handle, +void pcaFitTransform(const raft::handle_t& handle, double* input, double* trans_input, double* components, @@ -55,28 +55,28 @@ void pcaFitTransform(raft::handle_t& handle, double* noise_vars, const paramsPCA& prms, bool flip_signs_based_on_U); -void pcaInverseTransform(raft::handle_t& handle, +void pcaInverseTransform(const raft::handle_t& handle, float* trans_input, float* components, float* singular_vals, float* mu, float* input, const paramsPCA& prms); -void pcaInverseTransform(raft::handle_t& handle, +void pcaInverseTransform(const raft::handle_t& handle, double* trans_input, double* components, double* singular_vals, double* mu, double* input, const paramsPCA& prms); -void pcaTransform(raft::handle_t& handle, +void pcaTransform(const raft::handle_t& handle, float* input, float* components, float* trans_input, float* singular_vals, float* mu, const paramsPCA& prms); -void pcaTransform(raft::handle_t& handle, +void pcaTransform(const raft::handle_t& handle, double* input, double* components, double* trans_input, diff --git a/cpp/include/cuml/decomposition/tsvd.hpp b/cpp/include/cuml/decomposition/tsvd.hpp index abbc47eb8a..333886d39c 100644 --- a/cpp/include/cuml/decomposition/tsvd.hpp +++ b/cpp/include/cuml/decomposition/tsvd.hpp @@ -1,5 +1,5 @@ /* - * SPDX-FileCopyrightText: Copyright (c) 2018-2025, NVIDIA CORPORATION. + * SPDX-FileCopyrightText: Copyright (c) 2018-2026, NVIDIA CORPORATION. * SPDX-License-Identifier: Apache-2.0 */ @@ -13,39 +13,39 @@ class handle_t; namespace ML { -void tsvdFit(raft::handle_t& handle, +void tsvdFit(const raft::handle_t& handle, float* input, float* components, float* singular_vals, const paramsTSVD& prms, bool flip_signs_based_on_U); -void tsvdFit(raft::handle_t& handle, +void tsvdFit(const raft::handle_t& handle, double* input, double* components, double* singular_vals, const paramsTSVD& prms, bool flip_signs_based_on_U); -void tsvdInverseTransform(raft::handle_t& handle, +void tsvdInverseTransform(const raft::handle_t& handle, float* trans_input, float* components, float* input, const paramsTSVD& prms); -void tsvdInverseTransform(raft::handle_t& handle, +void tsvdInverseTransform(const raft::handle_t& handle, double* trans_input, double* components, double* input, const paramsTSVD& prms); -void tsvdTransform(raft::handle_t& handle, +void tsvdTransform(const raft::handle_t& handle, float* input, float* components, float* trans_input, const paramsTSVD& prms); -void tsvdTransform(raft::handle_t& handle, +void tsvdTransform(const raft::handle_t& handle, double* input, double* components, double* trans_input, const paramsTSVD& prms); -void tsvdFitTransform(raft::handle_t& handle, +void tsvdFitTransform(const raft::handle_t& handle, float* input, float* trans_input, float* components, @@ -54,7 +54,7 @@ void tsvdFitTransform(raft::handle_t& handle, float* singular_vals, const paramsTSVD& prms, bool flip_signs_based_on_U); -void tsvdFitTransform(raft::handle_t& handle, +void tsvdFitTransform(const raft::handle_t& handle, double* input, double* trans_input, double* components, diff --git a/cpp/src/pca/pca.cu b/cpp/src/pca/pca.cu index 8c93a961ad..f228f362fe 100644 --- a/cpp/src/pca/pca.cu +++ b/cpp/src/pca/pca.cu @@ -1,17 +1,122 @@ /* - * SPDX-FileCopyrightText: Copyright (c) 2018-2025, NVIDIA CORPORATION. + * SPDX-FileCopyrightText: Copyright (c) 2018-2026, NVIDIA CORPORATION. * SPDX-License-Identifier: Apache-2.0 */ -#include "pca.cuh" - #include +#include #include +#include namespace ML { -void pcaFit(raft::handle_t& handle, +template +void pca_fit_impl(const raft::handle_t& handle, + math_t* input, + math_t* components, + math_t* explained_var, + math_t* explained_var_ratio, + math_t* singular_vals, + math_t* mu, + math_t* noise_vars, + const paramsPCA& prms, + bool flip_signs_based_on_U) +{ + auto raft_prms = to_raft_params(prms); + raft::linalg::pca_fit( + handle, + raft_prms, + raft::make_device_matrix_view( + input, prms.n_rows, prms.n_cols), + raft::make_device_matrix_view( + components, prms.n_components, prms.n_cols), + raft::make_device_vector_view(explained_var, prms.n_components), + raft::make_device_vector_view(explained_var_ratio, prms.n_components), + raft::make_device_vector_view(singular_vals, prms.n_components), + raft::make_device_vector_view(mu, prms.n_cols), + raft::make_device_scalar_view(noise_vars), + flip_signs_based_on_U); +} + +template +void pca_fit_transform_impl(const raft::handle_t& handle, + math_t* input, + math_t* trans_input, + math_t* components, + math_t* explained_var, + math_t* explained_var_ratio, + math_t* singular_vals, + math_t* mu, + math_t* noise_vars, + const paramsPCA& prms, + bool flip_signs_based_on_U) +{ + auto raft_prms = to_raft_params(prms); + raft::linalg::pca_fit_transform( + handle, + raft_prms, + raft::make_device_matrix_view( + input, prms.n_rows, prms.n_cols), + raft::make_device_matrix_view( + trans_input, prms.n_rows, prms.n_components), + raft::make_device_matrix_view( + components, prms.n_components, prms.n_cols), + raft::make_device_vector_view(explained_var, prms.n_components), + raft::make_device_vector_view(explained_var_ratio, prms.n_components), + raft::make_device_vector_view(singular_vals, prms.n_components), + raft::make_device_vector_view(mu, prms.n_cols), + raft::make_device_scalar_view(noise_vars), + flip_signs_based_on_U); +} + +template +void pca_inverse_transform_impl(const raft::handle_t& handle, + math_t* trans_input, + math_t* components, + math_t* singular_vals, + math_t* mu, + math_t* input, + const paramsPCA& prms) +{ + auto raft_prms = to_raft_params(prms); + raft::linalg::pca_inverse_transform( + handle, + raft_prms, + raft::make_device_matrix_view( + trans_input, prms.n_rows, prms.n_components), + raft::make_device_matrix_view( + components, prms.n_components, prms.n_cols), + raft::make_device_vector_view(singular_vals, prms.n_components), + raft::make_device_vector_view(mu, prms.n_cols), + raft::make_device_matrix_view( + input, prms.n_rows, prms.n_cols)); +} + +template +void pca_transform_impl(const raft::handle_t& handle, + math_t* input, + math_t* components, + math_t* trans_input, + math_t* singular_vals, + math_t* mu, + const paramsPCA& prms) +{ + auto raft_prms = to_raft_params(prms); + raft::linalg::pca_transform( + handle, + raft_prms, + raft::make_device_matrix_view( + input, prms.n_rows, prms.n_cols), + raft::make_device_matrix_view( + components, prms.n_components, prms.n_cols), + raft::make_device_vector_view(singular_vals, prms.n_components), + raft::make_device_vector_view(mu, prms.n_cols), + raft::make_device_matrix_view( + trans_input, prms.n_rows, prms.n_components)); +} + +void pcaFit(const raft::handle_t& handle, float* input, float* components, float* explained_var, @@ -20,22 +125,21 @@ void pcaFit(raft::handle_t& handle, float* mu, float* noise_vars, const paramsPCA& prms, - bool flip_signs_based_on_U = false) + bool flip_signs_based_on_U) { - pcaFit(handle, - input, - components, - explained_var, - explained_var_ratio, - singular_vals, - mu, - noise_vars, - prms, - handle.get_stream(), - flip_signs_based_on_U); + pca_fit_impl(handle, + input, + components, + explained_var, + explained_var_ratio, + singular_vals, + mu, + noise_vars, + prms, + flip_signs_based_on_U); } -void pcaFit(raft::handle_t& handle, +void pcaFit(const raft::handle_t& handle, double* input, double* components, double* explained_var, @@ -44,22 +148,21 @@ void pcaFit(raft::handle_t& handle, double* mu, double* noise_vars, const paramsPCA& prms, - bool flip_signs_based_on_U = false) + bool flip_signs_based_on_U) { - pcaFit(handle, - input, - components, - explained_var, - explained_var_ratio, - singular_vals, - mu, - noise_vars, - prms, - handle.get_stream(), - flip_signs_based_on_U); + pca_fit_impl(handle, + input, + components, + explained_var, + explained_var_ratio, + singular_vals, + mu, + noise_vars, + prms, + flip_signs_based_on_U); } -void pcaFitTransform(raft::handle_t& handle, +void pcaFitTransform(const raft::handle_t& handle, float* input, float* trans_input, float* components, @@ -69,23 +172,22 @@ void pcaFitTransform(raft::handle_t& handle, float* mu, float* noise_vars, const paramsPCA& prms, - bool flip_signs_based_on_U = false) + bool flip_signs_based_on_U) { - pcaFitTransform(handle, - input, - trans_input, - components, - explained_var, - explained_var_ratio, - singular_vals, - mu, - noise_vars, - prms, - handle.get_stream(), - flip_signs_based_on_U); + pca_fit_transform_impl(handle, + input, + trans_input, + components, + explained_var, + explained_var_ratio, + singular_vals, + mu, + noise_vars, + prms, + flip_signs_based_on_U); } -void pcaFitTransform(raft::handle_t& handle, +void pcaFitTransform(const raft::handle_t& handle, double* input, double* trans_input, double* components, @@ -95,23 +197,22 @@ void pcaFitTransform(raft::handle_t& handle, double* mu, double* noise_vars, const paramsPCA& prms, - bool flip_signs_based_on_U = false) + bool flip_signs_based_on_U) { - pcaFitTransform(handle, - input, - trans_input, - components, - explained_var, - explained_var_ratio, - singular_vals, - mu, - noise_vars, - prms, - handle.get_stream(), - flip_signs_based_on_U); + pca_fit_transform_impl(handle, + input, + trans_input, + components, + explained_var, + explained_var_ratio, + singular_vals, + mu, + noise_vars, + prms, + flip_signs_based_on_U); } -void pcaInverseTransform(raft::handle_t& handle, +void pcaInverseTransform(const raft::handle_t& handle, float* trans_input, float* components, float* singular_vals, @@ -119,11 +220,10 @@ void pcaInverseTransform(raft::handle_t& handle, float* input, const paramsPCA& prms) { - pcaInverseTransform( - handle, trans_input, components, singular_vals, mu, input, prms, handle.get_stream()); + pca_inverse_transform_impl(handle, trans_input, components, singular_vals, mu, input, prms); } -void pcaInverseTransform(raft::handle_t& handle, +void pcaInverseTransform(const raft::handle_t& handle, double* trans_input, double* components, double* singular_vals, @@ -131,11 +231,10 @@ void pcaInverseTransform(raft::handle_t& handle, double* input, const paramsPCA& prms) { - pcaInverseTransform( - handle, trans_input, components, singular_vals, mu, input, prms, handle.get_stream()); + pca_inverse_transform_impl(handle, trans_input, components, singular_vals, mu, input, prms); } -void pcaTransform(raft::handle_t& handle, +void pcaTransform(const raft::handle_t& handle, float* input, float* components, float* trans_input, @@ -143,11 +242,10 @@ void pcaTransform(raft::handle_t& handle, float* mu, const paramsPCA& prms) { - pcaTransform( - handle, input, components, trans_input, singular_vals, mu, prms, handle.get_stream()); + pca_transform_impl(handle, input, components, trans_input, singular_vals, mu, prms); } -void pcaTransform(raft::handle_t& handle, +void pcaTransform(const raft::handle_t& handle, double* input, double* components, double* trans_input, @@ -155,8 +253,7 @@ void pcaTransform(raft::handle_t& handle, double* mu, const paramsPCA& prms) { - pcaTransform( - handle, input, components, trans_input, singular_vals, mu, prms, handle.get_stream()); + pca_transform_impl(handle, input, components, trans_input, singular_vals, mu, prms); } }; // end namespace ML diff --git a/cpp/src/pca/pca.cuh b/cpp/src/pca/pca.cuh deleted file mode 100644 index 1c713fa5fe..0000000000 --- a/cpp/src/pca/pca.cuh +++ /dev/null @@ -1,335 +0,0 @@ -/* - * SPDX-FileCopyrightText: Copyright (c) 2018-2026, NVIDIA CORPORATION. - * SPDX-License-Identifier: Apache-2.0 - */ - -#pragma once - -#include - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include - -#include - -namespace ML { - -template -void truncCompExpVars(const raft::handle_t& handle, - math_t* in, - math_t* components, - math_t* explained_var, - math_t* explained_var_ratio, - math_t* noise_vars, - const paramsTSVDTemplate& prms, - cudaStream_t stream) -{ - auto len = prms.n_cols * prms.n_cols; - rmm::device_uvector components_all(len, stream); - rmm::device_uvector explained_var_all(prms.n_cols, stream); - rmm::device_uvector explained_var_ratio_all(prms.n_cols, stream); - - calEig( - handle, in, components_all.data(), explained_var_all.data(), prms, stream); - raft::matrix::trunc_zero_origin( - handle, - raft::make_device_matrix_view( - components_all.data(), prms.n_cols, prms.n_cols), - raft::make_device_matrix_view( - components, prms.n_components, prms.n_cols)); - raft::matrix::ratio(handle, - raft::make_device_matrix_view( - explained_var_all.data(), prms.n_cols, std::size_t(1)), - raft::make_device_matrix_view( - explained_var_ratio_all.data(), prms.n_cols, std::size_t(1))); - raft::matrix::trunc_zero_origin( - handle, - raft::make_device_matrix_view( - explained_var_all.data(), prms.n_cols, std::size_t(1)), - raft::make_device_matrix_view( - explained_var, prms.n_components, std::size_t(1))); - raft::matrix::trunc_zero_origin( - handle, - raft::make_device_matrix_view( - explained_var_ratio_all.data(), prms.n_cols, std::size_t(1)), - raft::make_device_matrix_view( - explained_var_ratio, prms.n_components, std::size_t(1))); - - // Compute the scalar noise_vars defined as (pseudocode) - // (n_components < min(n_cols, n_rows)) ? explained_var_all[n_components:].mean() : 0 - if (prms.n_components < prms.n_cols && prms.n_components < prms.n_rows) { - raft::stats::mean(noise_vars, - explained_var_all.data() + prms.n_components, - std::size_t{1}, - prms.n_cols - prms.n_components, - false, - stream); - } else { - raft::matrix::fill( - handle, - raft::make_device_vector_view(noise_vars, std::size_t(1)), - math_t{0}); - } -} - -/** - * @brief perform fit operation for the pca. Generates eigenvectors, explained vars, singular vals, - * etc. - * @param[in] handle: cuml handle object - * @param[in] input: the data is fitted to PCA. Size n_rows x n_cols. The size of the data is - * indicated in prms. - * @param[out] components: the principal components of the input data. Size n_cols * n_components. - * @param[out] explained_var: explained variances (eigenvalues) of the principal components. Size - * n_components * 1. - * @param[out] explained_var_ratio: the ratio of the explained variance and total variance. Size - * n_components * 1. - * @param[out] singular_vals: singular values of the data. Size n_components * 1 - * @param[out] mu: mean of all the features (all the columns in the data). Size n_cols * 1. - * @param[out] noise_vars: variance of the noise. Size 1 * 1 (scalar). - * @param[in] prms: data structure that includes all the parameters from input size to algorithm. - * @param[in] stream cuda stream - */ -template -void pcaFit(const raft::handle_t& handle, - math_t* input, - math_t* components, - math_t* explained_var, - math_t* explained_var_ratio, - math_t* singular_vals, - math_t* mu, - math_t* noise_vars, - const paramsPCA& prms, - cudaStream_t stream, - bool flip_signs_based_on_U = false) -{ - auto cublas_handle = handle.get_cublas_handle(); - - ASSERT(prms.n_cols > 1, "Parameter n_cols: number of columns cannot be less than two"); - ASSERT(prms.n_rows > 1, "Parameter n_rows: number of rows cannot be less than two"); - ASSERT(prms.n_components > 0, - "Parameter n_components: number of components cannot be less than one"); - - auto n_components = prms.n_components; - if (n_components > prms.n_cols) n_components = prms.n_cols; - - raft::stats::mean(mu, input, prms.n_cols, prms.n_rows, false, stream); - - auto len = prms.n_cols * prms.n_cols; - rmm::device_uvector cov(len, stream); - - raft::stats::cov( - handle, cov.data(), input, mu, prms.n_cols, prms.n_rows, true, true, stream); - truncCompExpVars( - handle, cov.data(), components, explained_var, explained_var_ratio, noise_vars, prms, stream); - - math_t scalar = (prms.n_rows - 1); - raft::matrix::weighted_sqrt( - handle, - raft::make_device_matrix_view( - explained_var, std::size_t(1), n_components), - raft::make_device_matrix_view( - singular_vals, std::size_t(1), n_components), - raft::make_host_scalar_view(&scalar), - true); - - raft::stats::meanAdd(input, input, mu, prms.n_cols, prms.n_rows, stream); - - signFlipComponents(handle, - input, - components, - prms.n_rows, - prms.n_cols, - prms.n_components, - stream, - true, - flip_signs_based_on_U); -} - -/** - * @brief perform fit and transform operations for the pca. Generates transformed data, - * eigenvectors, explained vars, singular vals, etc. - * @param[in] handle: cuml handle object - * @param[in] input: the data is fitted to PCA. Size n_rows x n_cols. The size of the data is - * indicated in prms. - * @param[out] trans_input: the transformed data. Size n_rows * n_components. - * @param[out] components: the principal components of the input data. Size n_cols * n_components. - * @param[out] explained_var: explained variances (eigenvalues) of the principal components. Size - * n_components * 1. - * @param[out] explained_var_ratio: the ratio of the explained variance and total variance. Size - * n_components * 1. - * @param[out] singular_vals: singular values of the data. Size n_components * 1 - * @param[out] mu: mean of all the features (all the columns in the data). Size n_cols * 1. - * @param[out] noise_vars: variance of the noise. Size 1 * 1 (scalar). - * @param[in] prms: data structure that includes all the parameters from input size to algorithm. - * @param[in] stream cuda stream - */ -template -void pcaFitTransform(const raft::handle_t& handle, - math_t* input, - math_t* trans_input, - math_t* components, - math_t* explained_var, - math_t* explained_var_ratio, - math_t* singular_vals, - math_t* mu, - math_t* noise_vars, - const paramsPCA& prms, - cudaStream_t stream, - bool flip_signs_based_on_U = false) -{ - pcaFit(handle, - input, - components, - explained_var, - explained_var_ratio, - singular_vals, - mu, - noise_vars, - prms, - stream, - flip_signs_based_on_U); - pcaTransform(handle, input, components, trans_input, singular_vals, mu, prms, stream); -} - -// TODO: implement pcaGetCovariance function -template -void pcaGetCovariance() -{ - ASSERT(false, "pcaGetCovariance: will be implemented!"); -} - -// TODO: implement pcaGetPrecision function -template -void pcaGetPrecision() -{ - ASSERT(false, "pcaGetPrecision: will be implemented!"); -} - -/** - * @brief performs inverse transform operation for the pca. Transforms the transformed data back to - * original data. - * @param[in] handle: the internal cuml handle object - * @param[in] trans_input: the data is fitted to PCA. Size n_rows x n_components. - * @param[in] components: transpose of the principal components of the input data. Size n_components - * * n_cols. - * @param[in] singular_vals: singular values of the data. Size n_components * 1 - * @param[in] mu: mean of features (every column). - * @param[out] input: the data is fitted to PCA. Size n_rows x n_cols. - * @param[in] prms: data structure that includes all the parameters from input size to algorithm. - * @param[in] stream cuda stream - */ -template -void pcaInverseTransform(const raft::handle_t& handle, - math_t* trans_input, - math_t* components, - math_t* singular_vals, - math_t* mu, - math_t* input, - const paramsPCA& prms, - cudaStream_t stream) -{ - ASSERT(prms.n_cols > 1, "Parameter n_cols: number of columns cannot be less than two"); - ASSERT(prms.n_rows > 0, "Parameter n_rows: number of rows cannot be less than one"); - ASSERT(prms.n_components > 0, - "Parameter n_components: number of components cannot be less than one"); - - auto components_len = prms.n_cols * prms.n_components; - rmm::device_uvector components_copy{components_len, stream}; - raft::copy(components_copy.data(), components, prms.n_cols * prms.n_components, stream); - - if (prms.whiten) { - math_t sqrt_n_samples = sqrt(prms.n_rows - 1); - math_t scalar = prms.n_rows - 1 > 0 ? math_t(1 / sqrt_n_samples) : 0; - raft::linalg::scalarMultiply(components_copy.data(), - components_copy.data(), - scalar, - prms.n_cols * prms.n_components, - stream); - raft::linalg::binary_mult_skip_zero( - handle, - raft::make_device_matrix_view( - components_copy.data(), prms.n_cols, prms.n_components), - raft::make_device_vector_view(singular_vals, prms.n_components)); - } - - tsvdInverseTransform(handle, trans_input, components_copy.data(), input, prms, stream); - raft::stats::meanAdd(input, input, mu, prms.n_cols, prms.n_rows, stream); -} - -// TODO: implement pcaScore function -template -void pcaScore() -{ - ASSERT(false, "pcaScore: will be implemented!"); -} - -// TODO: implement pcaScoreSamples function -template -void pcaScoreSamples() -{ - ASSERT(false, "pcaScoreSamples: will be implemented!"); -} - -/** - * @brief performs transform operation for the pca. Transforms the data to eigenspace. - * @param[in] handle: the internal cuml handle object - * @param[in] input: the data is transformed. Size n_rows x n_components. - * @param[in] components: principal components of the input data. Size n_cols * n_components. - * @param[out] trans_input: the transformed data. Size n_rows * n_components. - * @param[in] singular_vals: singular values of the data. Size n_components * 1. - * @param[in] mu: mean value of the input data - * @param[in] prms: data structure that includes all the parameters from input size to algorithm. - * @param[in] stream cuda stream - */ -template -void pcaTransform(const raft::handle_t& handle, - math_t* input, - math_t* components, - math_t* trans_input, - math_t* singular_vals, - math_t* mu, - const paramsPCA& prms, - cudaStream_t stream) -{ - ASSERT(prms.n_cols > 1, "Parameter n_cols: number of columns cannot be less than two"); - ASSERT(prms.n_rows > 0, "Parameter n_rows: number of rows cannot be less than one"); - ASSERT(prms.n_components > 0, - "Parameter n_components: number of components cannot be less than one"); - - auto components_len = prms.n_cols * prms.n_components; - rmm::device_uvector components_copy{components_len, stream}; - raft::copy(components_copy.data(), components, prms.n_cols * prms.n_components, stream); - - if (prms.whiten) { - math_t scalar = math_t(sqrt(prms.n_rows - 1)); - raft::linalg::scalarMultiply(components_copy.data(), - components_copy.data(), - scalar, - prms.n_cols * prms.n_components, - stream); - raft::linalg::binary_div_skip_zero( - handle, - raft::make_device_matrix_view( - components_copy.data(), prms.n_cols, prms.n_components), - raft::make_device_vector_view(singular_vals, prms.n_components)); - } - - raft::stats::meanCenter(input, input, mu, prms.n_cols, prms.n_rows, stream); - tsvdTransform(handle, input, components_copy.data(), trans_input, prms, stream); - raft::stats::meanAdd(input, input, mu, prms.n_cols, prms.n_rows, stream); -} - -}; // end namespace ML diff --git a/cpp/src/pca/pca_mg.cu b/cpp/src/pca/pca_mg.cu index 147c698938..325940c87e 100644 --- a/cpp/src/pca/pca_mg.cu +++ b/cpp/src/pca/pca_mg.cu @@ -3,8 +3,6 @@ * SPDX-License-Identifier: Apache-2.0 */ -#include "pca.cuh" - #include #include #include @@ -15,8 +13,12 @@ #include #include +#include +#include #include +#include #include +#include #include #include #include @@ -56,12 +58,23 @@ void fit_impl(raft::handle_t& handle, MLCommon::Stats::opg::cov(handle, cov, input_data, input_desc, mu_data, true, streams, n_streams); - ML::truncCompExpVars( - handle, cov.ptr, components, explained_var, explained_var_ratio, noise_vars, prms, streams[0]); - - T scalar = (prms.n_rows - 1); + auto raft_prms = to_raft_params(prms); raft::resources handle_stream_zero; raft::resource::set_cuda_stream(handle_stream_zero, streams[0]); + + raft::linalg::trunc_comp_exp_vars( + handle_stream_zero, + raft_prms, + raft::make_device_matrix_view( + cov.ptr, prms.n_cols, prms.n_cols), + raft::make_device_matrix_view( + components, prms.n_components, prms.n_cols), + raft::make_device_vector_view(explained_var, prms.n_components), + raft::make_device_vector_view(explained_var_ratio, prms.n_components), + raft::make_device_scalar_view(noise_vars), + prms.n_rows); + + T scalar = (prms.n_rows - 1); raft::matrix::weighted_sqrt(handle_stream_zero, raft::make_device_matrix_view( explained_var, std::size_t(1), prms.n_components), @@ -87,15 +100,14 @@ void fit_impl(raft::handle_t& handle, for (std::uint32_t i = 0; i < n_streams; i++) { handle.sync_stream(streams[i]); } - signFlipComponents(handle, - input_data[0]->ptr, - components, - prms.n_rows, - prms.n_cols, - prms.n_components, - streams[0], - true, - false); + raft::linalg::sign_flip_components( + handle_stream_zero, + raft::make_device_matrix_view( + input_data[0]->ptr, prms.n_rows, prms.n_cols), + raft::make_device_matrix_view( + components, prms.n_components, prms.n_cols), + true, + false); } } diff --git a/cpp/src/tsne/tsne_runner.cuh b/cpp/src/tsne/tsne_runner.cuh index a54e56e7fa..9493735e92 100644 --- a/cpp/src/tsne/tsne_runner.cuh +++ b/cpp/src/tsne/tsne_runner.cuh @@ -1,5 +1,5 @@ /* - * SPDX-FileCopyrightText: Copyright (c) 2019-2025, NVIDIA CORPORATION. + * SPDX-FileCopyrightText: Copyright (c) 2019-2026, NVIDIA CORPORATION. * SPDX-License-Identifier: Apache-2.0 */ @@ -13,12 +13,15 @@ #include "utils.cuh" #include +#include #include #include #include #include #include +#include +#include #include #include @@ -26,7 +29,6 @@ #include #include -#include #include #include @@ -116,7 +118,7 @@ class TSNE_runner { mu.data_handle(), noise_vars.data_handle(), prms, - stream); + false); auto mean_result = raft::make_device_vector(handle, dim); auto stddev_result = raft::make_device_vector(handle, dim); diff --git a/cpp/src/tsvd/tsvd.cu b/cpp/src/tsvd/tsvd.cu index d5ae22f2de..56049cda02 100644 --- a/cpp/src/tsvd/tsvd.cu +++ b/cpp/src/tsvd/tsvd.cu @@ -1,39 +1,125 @@ /* - * SPDX-FileCopyrightText: Copyright (c) 2018-2025, NVIDIA CORPORATION. + * SPDX-FileCopyrightText: Copyright (c) 2018-2026, NVIDIA CORPORATION. * SPDX-License-Identifier: Apache-2.0 */ -#include "tsvd.cuh" - #include +#include #include +#include namespace ML { -void tsvdFit(raft::handle_t& handle, +namespace { + +template +void tsvd_fit_impl(const raft::handle_t& handle, + math_t* input, + math_t* components, + math_t* singular_vals, + const paramsTSVD& prms, + bool flip_signs_based_on_U) +{ + auto raft_prms = to_raft_params(prms); + raft::linalg::tsvd_fit( + handle, + raft_prms, + raft::make_device_matrix_view( + input, prms.n_rows, prms.n_cols), + raft::make_device_matrix_view( + components, prms.n_components, prms.n_cols), + raft::make_device_vector_view(singular_vals, prms.n_components), + flip_signs_based_on_U); +} + +template +void tsvd_fit_transform_impl(const raft::handle_t& handle, + math_t* input, + math_t* trans_input, + math_t* components, + math_t* explained_var, + math_t* explained_var_ratio, + math_t* singular_vals, + const paramsTSVD& prms, + bool flip_signs_based_on_U) +{ + auto raft_prms = to_raft_params(prms); + raft::linalg::tsvd_fit_transform( + handle, + raft_prms, + raft::make_device_matrix_view( + input, prms.n_rows, prms.n_cols), + raft::make_device_matrix_view( + trans_input, prms.n_rows, prms.n_components), + raft::make_device_matrix_view( + components, prms.n_components, prms.n_cols), + raft::make_device_vector_view(explained_var, prms.n_components), + raft::make_device_vector_view(explained_var_ratio, prms.n_components), + raft::make_device_vector_view(singular_vals, prms.n_components), + flip_signs_based_on_U); +} + +template +void tsvd_transform_impl(const raft::handle_t& handle, + math_t* input, + math_t* components, + math_t* trans_input, + const paramsTSVD& prms) +{ + auto raft_prms = to_raft_params(prms); + raft::linalg::tsvd_transform(handle, + raft_prms, + raft::make_device_matrix_view( + input, prms.n_rows, prms.n_cols), + raft::make_device_matrix_view( + components, prms.n_components, prms.n_cols), + raft::make_device_matrix_view( + trans_input, prms.n_rows, prms.n_components)); +} + +template +void tsvd_inverse_transform_impl(const raft::handle_t& handle, + math_t* trans_input, + math_t* components, + math_t* input, + const paramsTSVD& prms) +{ + auto raft_prms = to_raft_params(prms); + raft::linalg::tsvd_inverse_transform( + handle, + raft_prms, + raft::make_device_matrix_view( + trans_input, prms.n_rows, prms.n_components), + raft::make_device_matrix_view( + components, prms.n_components, prms.n_cols), + raft::make_device_matrix_view( + input, prms.n_rows, prms.n_cols)); +} + +} // anonymous namespace + +void tsvdFit(const raft::handle_t& handle, float* input, float* components, float* singular_vals, const paramsTSVD& prms, - bool flip_signs_based_on_U = false) + bool flip_signs_based_on_U) { - tsvdFit( - handle, input, components, singular_vals, prms, handle.get_stream(), flip_signs_based_on_U); + tsvd_fit_impl(handle, input, components, singular_vals, prms, flip_signs_based_on_U); } -void tsvdFit(raft::handle_t& handle, +void tsvdFit(const raft::handle_t& handle, double* input, double* components, double* singular_vals, const paramsTSVD& prms, - bool flip_signs_based_on_U = false) + bool flip_signs_based_on_U) { - tsvdFit( - handle, input, components, singular_vals, prms, handle.get_stream(), flip_signs_based_on_U); + tsvd_fit_impl(handle, input, components, singular_vals, prms, flip_signs_based_on_U); } -void tsvdFitTransform(raft::handle_t& handle, +void tsvdFitTransform(const raft::handle_t& handle, float* input, float* trans_input, float* components, @@ -41,21 +127,20 @@ void tsvdFitTransform(raft::handle_t& handle, float* explained_var_ratio, float* singular_vals, const paramsTSVD& prms, - bool flip_signs_based_on_U = false) + bool flip_signs_based_on_U) { - tsvdFitTransform(handle, - input, - trans_input, - components, - explained_var, - explained_var_ratio, - singular_vals, - prms, - handle.get_stream(), - flip_signs_based_on_U); + tsvd_fit_transform_impl(handle, + input, + trans_input, + components, + explained_var, + explained_var_ratio, + singular_vals, + prms, + flip_signs_based_on_U); } -void tsvdFitTransform(raft::handle_t& handle, +void tsvdFitTransform(const raft::handle_t& handle, double* input, double* trans_input, double* components, @@ -63,54 +148,53 @@ void tsvdFitTransform(raft::handle_t& handle, double* explained_var_ratio, double* singular_vals, const paramsTSVD& prms, - bool flip_signs_based_on_U = false) + bool flip_signs_based_on_U) { - tsvdFitTransform(handle, - input, - trans_input, - components, - explained_var, - explained_var_ratio, - singular_vals, - prms, - handle.get_stream(), - flip_signs_based_on_U); + tsvd_fit_transform_impl(handle, + input, + trans_input, + components, + explained_var, + explained_var_ratio, + singular_vals, + prms, + flip_signs_based_on_U); } -void tsvdTransform(raft::handle_t& handle, +void tsvdTransform(const raft::handle_t& handle, float* input, float* components, float* trans_input, const paramsTSVD& prms) { - tsvdTransform(handle, input, components, trans_input, prms, handle.get_stream()); + tsvd_transform_impl(handle, input, components, trans_input, prms); } -void tsvdTransform(raft::handle_t& handle, +void tsvdTransform(const raft::handle_t& handle, double* input, double* components, double* trans_input, const paramsTSVD& prms) { - tsvdTransform(handle, input, components, trans_input, prms, handle.get_stream()); + tsvd_transform_impl(handle, input, components, trans_input, prms); } -void tsvdInverseTransform(raft::handle_t& handle, +void tsvdInverseTransform(const raft::handle_t& handle, float* trans_input, float* components, float* input, const paramsTSVD& prms) { - tsvdInverseTransform(handle, trans_input, components, input, prms, handle.get_stream()); + tsvd_inverse_transform_impl(handle, trans_input, components, input, prms); } -void tsvdInverseTransform(raft::handle_t& handle, +void tsvdInverseTransform(const raft::handle_t& handle, double* trans_input, double* components, double* input, const paramsTSVD& prms) { - tsvdInverseTransform(handle, trans_input, components, input, prms, handle.get_stream()); + tsvd_inverse_transform_impl(handle, trans_input, components, input, prms); } }; // end namespace ML diff --git a/cpp/src/tsvd/tsvd.cuh b/cpp/src/tsvd/tsvd.cuh deleted file mode 100644 index ce73567301..0000000000 --- a/cpp/src/tsvd/tsvd.cuh +++ /dev/null @@ -1,509 +0,0 @@ -/* - * SPDX-FileCopyrightText: Copyright (c) 2018-2026, NVIDIA CORPORATION. - * SPDX-License-Identifier: Apache-2.0 - */ - -#pragma once - -#include - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include -#include -#include - -#include -#include -#include -#include - -namespace ML { - -template -void calCompExpVarsSvd(const raft::handle_t& handle, - math_t* in, - math_t* components, - math_t* singular_vals, - math_t* explained_vars, - math_t* explained_var_ratio, - const paramsTSVD& prms, - cudaStream_t stream) -{ - auto cusolver_handle = handle.get_cusolver_dn_handle(); - auto cublas_handle = handle.get_cublas_handle(); - - auto diff = prms.n_cols - prms.n_components; - math_t ratio = math_t(diff) / math_t(prms.n_cols); - ASSERT(ratio >= math_t(0.2), - "Number of components should be less than at least 80 percent of the " - "number of features"); - - std::size_t p = static_cast(math_t(0.1) * math_t(prms.n_cols)); - // int p = int(math_t(prms.n_cols) / math_t(4)); - ASSERT(p >= 5, "RSVD should be used where the number of columns are at least 50"); - - auto total_random_vecs = prms.n_components + p; - ASSERT(total_random_vecs < prms.n_cols, - "RSVD should be used where the number of columns are at least 50"); - - rmm::device_uvector components_temp(prms.n_cols * prms.n_components, stream); - math_t* left_eigvec = nullptr; - raft::linalg::rsvdFixedRank(handle, - in, - prms.n_rows, - prms.n_cols, - singular_vals, - left_eigvec, - components_temp.data(), - prms.n_components, - p, - true, - false, - true, - false, - (math_t)prms.tol, - prms.n_iterations, - stream); - - raft::linalg::transpose( - handle, components_temp.data(), components, prms.n_cols, prms.n_components, stream); - - raft::matrix::weighted_power( - handle, - raft::make_device_matrix_view( - singular_vals, std::size_t(1), prms.n_components), - raft::make_device_matrix_view( - explained_vars, std::size_t(1), prms.n_components), - math_t(1)); - raft::matrix::ratio(handle, explained_vars, explained_var_ratio, prms.n_components, stream); -} - -template -void calEig(const raft::handle_t& handle, - math_t* in, - math_t* components, - math_t* explained_var, - const paramsTSVDTemplate& prms, - cudaStream_t stream) -{ - auto cusolver_handle = handle.get_cusolver_dn_handle(); - - if (prms.algorithm == enum_solver::COV_EIG_JACOBI) { - raft::linalg::eigJacobi(handle, - in, - prms.n_cols, - prms.n_cols, - components, - explained_var, - stream, - (math_t)prms.tol, - prms.n_iterations); - } else { - raft::linalg::eigDC(handle, in, prms.n_cols, prms.n_cols, components, explained_var, stream); - } - raft::resources handle_stream_zero; - raft::resource::set_cuda_stream(handle_stream_zero, stream); - - raft::matrix::col_reverse(handle_stream_zero, - raft::make_device_matrix_view( - components, prms.n_cols, prms.n_cols)); - raft::linalg::transpose(components, prms.n_cols, stream); - - raft::matrix::row_reverse(handle_stream_zero, - raft::make_device_matrix_view( - explained_var, prms.n_cols, std::size_t(1))); -} - -/** - * @defgroup sign flip for PCA and tSVD. This is used to stabilize the sign of column major eigen - * vectors - * @param handle: resource handle - * @param components: components matrix, used to determine the sign of max absolute value - * @param input: input data - * @param n_rows: number of rows of components matrix - * @param n_cols: number of columns of components matrix - * @param n_samples: number of samples (number of rows of input) - * @param stream: cuda stream - * @param flip_signs_based_on_U whether to determine signs by U (true) or V.T (false) - * @{ - */ -template -void signFlipComponents(const raft::handle_t& handle, - math_t* input, - math_t* components, - std::size_t n_samples, - std::size_t n_features, - std::size_t n_components, - cudaStream_t stream, - bool center, - bool flip_signs_based_on_U = false) -{ - rmm::device_uvector max_vals(n_components, stream); - auto components_view = raft::make_device_matrix_view( - components, n_components, n_features); - auto max_vals_view = - raft::make_device_vector_view(max_vals.data(), n_components); - - // Step 1: find U or V max absolute values - // X = U @ S @ V - // X: input matrix, n_samples * n_features - // U: n_samples * n_components - // S: diagonal matrix of eigen-values, n_components * n_components - // V: components, n_components * n_features - // U @ S = X @ V.T, where the signs of U @ S are solely determined by U - if (flip_signs_based_on_U) { - if (center) { - // If center, X -= X.mean(axis=0) - rmm::device_uvector col_means(n_features, stream); - raft::stats::mean(col_means.data(), input, n_features, n_samples, stream); - raft::stats::meanCenter( - input, input, col_means.data(), n_features, n_samples, stream); - } - rmm::device_uvector US(n_samples * n_components, stream); - raft::linalg::gemm(handle, - input, - n_samples, - n_features, - components, - US.data(), - n_samples, - n_components, - CUBLAS_OP_N, - CUBLAS_OP_T, - math_t(1), - math_t(0), - stream); - raft::linalg::reduce( - max_vals.data(), - US.data(), - n_components, - n_samples, - math_t(0), - stream, - false, - raft::identity_op(), - [] __device__(math_t a, math_t b) { - math_t abs_a = a >= math_t(0) ? a : -a; - math_t abs_b = b >= math_t(0) ? b : -b; - return abs_a >= abs_b ? a : b; - }, - raft::identity_op()); - } else { - raft::linalg::reduce( - max_vals.data(), - components, - n_features, - n_components, - math_t(0), - stream, - false, - raft::identity_op(), - [] __device__(math_t a, math_t b) { - math_t abs_a = a >= math_t(0) ? a : -a; - math_t abs_b = b >= math_t(0) ? b : -b; - return abs_a >= abs_b ? a : b; - }, - raft::identity_op()); - } - - // Step 2: flip rows where needed - raft::linalg::map_offset( - handle, - components_view, - [components_view, max_vals_view, n_components, n_features] __device__(auto idx) { - std::size_t row = idx % n_components; - std::size_t column = idx / n_components; - return (max_vals_view(row) < math_t(0)) ? (-components_view(row, column)) - : components_view(row, column); - }); -} - -/** - * @defgroup sign flip for PCA and tSVD. This is used to stabilize the sign of column major eigen - * vectors - * @param input: input matrix that will be used to determine the sign. - * @param n_rows: number of rows of input matrix - * @param n_cols: number of columns of input matrix - * @param components: components matrix. - * @param n_cols_comp: number of columns of components matrix - * @param stream cuda stream - * @{ - */ -template -void signFlip(math_t* input, - std::size_t n_rows, - std::size_t n_cols, - math_t* components, - std::size_t n_cols_comp, - cudaStream_t stream) -{ - auto counting = thrust::make_counting_iterator(0); - auto m = n_rows; - - thrust::for_each( - rmm::exec_policy(stream), counting, counting + n_cols, [=] __device__(std::size_t idx) { - auto d_i = idx * m; - auto end = d_i + m; - - math_t max = 0.0; - std::size_t max_index = 0; - for (auto i = d_i; i < end; i++) { - math_t val = input[i]; - if (val < 0.0) { val = -val; } - if (val > max) { - max = val; - max_index = i; - } - } - - if (input[max_index] < 0.0) { - for (auto i = d_i; i < end; i++) { - input[i] = -input[i]; - } - - auto len = n_cols * n_cols_comp; - for (auto i = idx; i < len; i = i + n_cols) { - components[i] = -components[i]; - } - } - }); -} - -/** - * @brief perform fit operation for the tsvd. Generates eigenvectors, explained vars, singular vals, - * etc. - * @param[in] handle: the internal cuml handle object - * @param[in] input: the data is fitted to PCA. Size n_rows x n_cols. The size of the data is - * indicated in prms. - * @param[out] components: the principal components of the input data. Size n_cols * n_components. - * @param[out] singular_vals: singular values of the data. Size n_components * 1 - * @param[in] prms: data structure that includes all the parameters from input size to algorithm. - * @param[in] stream cuda stream - */ -template -void tsvdFit(const raft::handle_t& handle, - math_t* input, - math_t* components, - math_t* singular_vals, - const paramsTSVD& prms, - cudaStream_t stream, - bool flip_signs_based_on_U = false) -{ - auto cublas_handle = handle.get_cublas_handle(); - - ASSERT(prms.n_cols > 1, "Parameter n_cols: number of columns cannot be less than two"); - ASSERT(prms.n_rows > 1, "Parameter n_rows: number of rows cannot be less than two"); - ASSERT(prms.n_components > 0, - "Parameter n_components: number of components cannot be less than one"); - - auto n_components = prms.n_components; - if (prms.n_components > prms.n_cols) n_components = prms.n_cols; - - size_t len = prms.n_cols * prms.n_cols; - rmm::device_uvector input_cross_mult(len, stream); - - math_t alpha = math_t(1); - math_t beta = math_t(0); - raft::linalg::gemm(handle, - input, - prms.n_rows, - prms.n_cols, - input, - input_cross_mult.data(), - prms.n_cols, - prms.n_cols, - CUBLAS_OP_T, - CUBLAS_OP_N, - alpha, - beta, - stream); - - rmm::device_uvector components_all(len, stream); - rmm::device_uvector explained_var_all(prms.n_cols, stream); - - calEig( - handle, input_cross_mult.data(), components_all.data(), explained_var_all.data(), prms, stream); - - raft::matrix::trunc_zero_origin( - handle, - raft::make_device_matrix_view( - components_all.data(), prms.n_cols, prms.n_cols), - raft::make_device_matrix_view( - components, n_components, prms.n_cols)); - - math_t scalar = math_t(1); - raft::matrix::weighted_sqrt( - handle, - raft::make_device_matrix_view( - explained_var_all.data(), std::size_t(1), n_components), - raft::make_device_matrix_view( - singular_vals, std::size_t(1), n_components), - raft::make_host_scalar_view(&scalar)); - - signFlipComponents(handle, - input, - components, - prms.n_rows, - prms.n_cols, - n_components, - stream, - false, - flip_signs_based_on_U); -} - -/** - * @brief performs fit and transform operations for the tsvd. Generates transformed data, - * eigenvectors, explained vars, singular vals, etc. - * @param[in] handle: the internal cuml handle object - * @param[in] input: the data is fitted to PCA. Size n_rows x n_cols. The size of the data is - * indicated in prms. - * @param[out] trans_input: the transformed data. Size n_rows * n_components. - * @param[out] components: the principal components of the input data. Size n_cols * n_components. - * @param[out] explained_var: explained variances (eigenvalues) of the principal components. Size - * n_components * 1. - * @param[out] explained_var_ratio: the ratio of the explained variance and total variance. Size - * n_components * 1. - * @param[out] singular_vals: singular values of the data. Size n_components * 1 - * @param[in] prms: data structure that includes all the parameters from input size to algorithm. - * @param[in] stream cuda stream - */ -template -void tsvdFitTransform(const raft::handle_t& handle, - math_t* input, - math_t* trans_input, - math_t* components, - math_t* explained_var, - math_t* explained_var_ratio, - math_t* singular_vals, - const paramsTSVD& prms, - cudaStream_t stream, - bool flip_signs_based_on_U = false) -{ - tsvdFit(handle, input, components, singular_vals, prms, stream, flip_signs_based_on_U); - tsvdTransform(handle, input, components, trans_input, prms, stream); - - rmm::device_uvector mu_trans(prms.n_components, stream); - raft::stats::mean( - mu_trans.data(), trans_input, prms.n_components, prms.n_rows, false, stream); - raft::stats::vars( - explained_var, trans_input, mu_trans.data(), prms.n_components, prms.n_rows, false, stream); - - rmm::device_uvector mu(prms.n_cols, stream); - rmm::device_uvector vars(prms.n_cols, stream); - - raft::stats::mean(mu.data(), input, prms.n_cols, prms.n_rows, false, stream); - raft::stats::vars(vars.data(), input, mu.data(), prms.n_cols, prms.n_rows, false, stream); - - rmm::device_scalar total_vars(stream); - raft::stats::sum(total_vars.data(), vars.data(), std::size_t(1), prms.n_cols, stream); - - math_t total_vars_h; - raft::update_host(&total_vars_h, total_vars.data(), 1, stream); - handle.sync_stream(stream); - math_t scalar = math_t(1) / total_vars_h; - - raft::linalg::scalarMultiply( - explained_var_ratio, explained_var, scalar, prms.n_components, stream); -} - -/** - * @brief performs transform operation for the tsvd. Transforms the data to eigenspace. - * @param[in] handle the internal cuml handle object - * @param[in] input: the data is transformed. Size n_rows x n_components. - * @param[in] components: principal components of the input data. Size n_cols * n_components. - * @param[out] trans_input: output that is transformed version of input - * @param[in] prms: data structure that includes all the parameters from input size to algorithm. - * @param[in] stream cuda stream - */ -template -void tsvdTransform(const raft::handle_t& handle, - math_t* input, - math_t* components, - math_t* trans_input, - const paramsTSVD& prms, - cudaStream_t stream) -{ - ASSERT(prms.n_cols > 1, "Parameter n_cols: number of columns cannot be less than two"); - ASSERT(prms.n_rows > 0, "Parameter n_rows: number of rows cannot be less than one"); - ASSERT(prms.n_components > 0, - "Parameter n_components: number of components cannot be less than one"); - - math_t alpha = math_t(1); - math_t beta = math_t(0); - raft::linalg::gemm(handle, - input, - prms.n_rows, - prms.n_cols, - components, - trans_input, - prms.n_rows, - prms.n_components, - CUBLAS_OP_N, - CUBLAS_OP_T, - alpha, - beta, - stream); -} - -/** - * @brief performs inverse transform operation for the tsvd. Transforms the transformed data back to - * original data. - * @param[in] handle the internal cuml handle object - * @param[in] trans_input: the data is fitted to PCA. Size n_rows x n_components. - * @param[in] components: transpose of the principal components of the input data. Size n_components - * * n_cols. - * @param[out] input: the data is fitted to PCA. Size n_rows x n_cols. - * @param[in] prms: data structure that includes all the parameters from input size to algorithm. - * @param[in] stream cuda stream - */ -template -void tsvdInverseTransform(const raft::handle_t& handle, - math_t* trans_input, - math_t* components, - math_t* input, - const paramsTSVD& prms, - cudaStream_t stream) -{ - ASSERT(prms.n_cols > 1, "Parameter n_cols: number of columns cannot be less than one"); - ASSERT(prms.n_rows > 0, "Parameter n_rows: number of rows cannot be less than one"); - ASSERT(prms.n_components > 0, - "Parameter n_components: number of components cannot be less than one"); - - math_t alpha = math_t(1); - math_t beta = math_t(0); - - raft::linalg::gemm(handle, - trans_input, - prms.n_rows, - prms.n_components, - components, - input, - prms.n_rows, - prms.n_cols, - CUBLAS_OP_N, - CUBLAS_OP_N, - alpha, - beta, - stream); -} - -}; // end namespace ML diff --git a/cpp/src/tsvd/tsvd_mg.cu b/cpp/src/tsvd/tsvd_mg.cu index dece19ecac..dcbe2cdcd1 100644 --- a/cpp/src/tsvd/tsvd_mg.cu +++ b/cpp/src/tsvd/tsvd_mg.cu @@ -3,8 +3,6 @@ * SPDX-License-Identifier: Apache-2.0 */ -#include "tsvd.cuh" - #include #include #include @@ -16,6 +14,7 @@ #include #include #include +#include #include #include #include @@ -52,10 +51,19 @@ void fit_impl(raft::handle_t& handle, rmm::device_uvector components_all(len, streams[0]); rmm::device_uvector explained_var_all(prms.n_cols, streams[0]); - ML::calEig(handle, cov.ptr, components_all.data(), explained_var_all.data(), prms, streams[0]); + auto raft_prms = to_raft_params(prms); raft::resources handle_stream_zero; raft::resource::set_cuda_stream(handle_stream_zero, streams[0]); + + raft::linalg::cal_eig( + handle_stream_zero, + raft_prms, + raft::make_device_matrix_view( + cov.ptr, prms.n_cols, prms.n_cols), + raft::make_device_matrix_view( + components_all.data(), prms.n_cols, prms.n_cols), + raft::make_device_vector_view(explained_var_all.data(), prms.n_cols)); raft::matrix::trunc_zero_origin( handle_stream_zero, raft::make_device_matrix_view( @@ -83,15 +91,14 @@ void fit_impl(raft::handle_t& handle, n_streams, false); } else { - signFlipComponents(handle, - input_data[0]->ptr, - components, - prms.n_rows, - prms.n_cols, - prms.n_components, - streams[0], - false, - false); + raft::linalg::sign_flip_components( + handle_stream_zero, + raft::make_device_matrix_view( + input_data[0]->ptr, prms.n_rows, prms.n_cols), + raft::make_device_matrix_view( + components, prms.n_components, prms.n_cols), + false, + false); } } diff --git a/cpp/tests/sg/pca_test.cu b/cpp/tests/sg/pca_test.cu index 1409bc85c9..fa12c9c983 100644 --- a/cpp/tests/sg/pca_test.cu +++ b/cpp/tests/sg/pca_test.cu @@ -1,9 +1,10 @@ /* - * SPDX-FileCopyrightText: Copyright (c) 2018-2025, NVIDIA CORPORATION. + * SPDX-FileCopyrightText: Copyright (c) 2018-2026, NVIDIA CORPORATION. * SPDX-License-Identifier: Apache-2.0 */ #include +#include #include #include @@ -11,7 +12,6 @@ #include #include -#include #include #include @@ -106,23 +106,21 @@ class PcaTest : public ::testing::TestWithParam> { mean.data(), noise_vars.data(), prms, - stream); + false); pcaTransform(handle, data.data(), components.data(), trans_data.data(), singular_vals.data(), mean.data(), - prms, - stream); + prms); pcaInverseTransform(handle, trans_data.data(), components.data(), singular_vals.data(), mean.data(), data_back.data(), - prms, - stream); + prms); } void advancedTest() @@ -161,7 +159,7 @@ class PcaTest : public ::testing::TestWithParam> { mean2.data(), noise_vars2.data(), prms, - stream); + false); pcaInverseTransform(handle, data2_trans.data(), @@ -169,8 +167,7 @@ class PcaTest : public ::testing::TestWithParam> { singular_vals2.data(), mean2.data(), data2_back.data(), - prms, - stream); + prms); } protected: diff --git a/cpp/tests/sg/tsvd_test.cu b/cpp/tests/sg/tsvd_test.cu index bacc837b79..2fdc39c46b 100644 --- a/cpp/tests/sg/tsvd_test.cu +++ b/cpp/tests/sg/tsvd_test.cu @@ -1,9 +1,10 @@ /* - * SPDX-FileCopyrightText: Copyright (c) 2018-2025, NVIDIA CORPORATION. + * SPDX-FileCopyrightText: Copyright (c) 2018-2026, NVIDIA CORPORATION. * SPDX-License-Identifier: Apache-2.0 */ #include +#include #include #include @@ -11,7 +12,6 @@ #include #include -#include #include @@ -82,7 +82,7 @@ class TsvdTest : public ::testing::TestWithParam> { else prms.algorithm = solver::COV_EIG_JACOBI; - tsvdFit(handle, data.data(), components.data(), singular_vals.data(), prms, stream); + tsvdFit(handle, data.data(), components.data(), singular_vals.data(), prms, false); } void advancedTest() @@ -130,11 +130,10 @@ class TsvdTest : public ::testing::TestWithParam> { explained_var_ratio2.data(), singular_vals2.data(), prms, - stream); + false); data2_back.resize(len, stream); - tsvdInverseTransform( - handle, data2_trans.data(), components2.data(), data2_back.data(), prms, stream); + tsvdInverseTransform(handle, data2_trans.data(), components2.data(), data2_back.data(), prms); } protected: