From afd395d9baa19f471f0c587d6c44c7b5e2b0bd3c Mon Sep 17 00:00:00 2001 From: aamijar Date: Fri, 13 Feb 2026 09:09:30 +0000 Subject: [PATCH 01/12] move-pca-from-cuml --- cpp/include/raft/linalg/detail/pca.cuh | 333 ++++++++++++++++ cpp/include/raft/linalg/detail/tsvd.cuh | 508 ++++++++++++++++++++++++ cpp/include/raft/linalg/pca.cuh | 194 +++++++++ cpp/include/raft/linalg/pca_types.hpp | 78 ++++ cpp/include/raft/linalg/tsvd.cuh | 199 ++++++++++ cpp/tests/CMakeLists.txt | 4 +- cpp/tests/linalg/pca.cu | 316 +++++++++++++++ cpp/tests/linalg/tsvd.cu | 206 ++++++++++ 8 files changed, 1837 insertions(+), 1 deletion(-) create mode 100644 cpp/include/raft/linalg/detail/pca.cuh create mode 100644 cpp/include/raft/linalg/detail/tsvd.cuh create mode 100644 cpp/include/raft/linalg/pca.cuh create mode 100644 cpp/include/raft/linalg/pca_types.hpp create mode 100644 cpp/include/raft/linalg/tsvd.cuh create mode 100644 cpp/tests/linalg/pca.cu create mode 100644 cpp/tests/linalg/tsvd.cu diff --git a/cpp/include/raft/linalg/detail/pca.cuh b/cpp/include/raft/linalg/detail/pca.cuh new file mode 100644 index 0000000000..43d7640efd --- /dev/null +++ b/cpp/include/raft/linalg/detail/pca.cuh @@ -0,0 +1,333 @@ +/* + * 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 raft::linalg::detail { + +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); + + detail::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); + detail::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); + + detail::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) +{ + detail::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)); + } + + detail::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); + detail::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 raft::linalg::detail diff --git a/cpp/include/raft/linalg/detail/tsvd.cuh b/cpp/include/raft/linalg/detail/tsvd.cuh new file mode 100644 index 0000000000..39adffc8e3 --- /dev/null +++ b/cpp/include/raft/linalg/detail/tsvd.cuh @@ -0,0 +1,508 @@ +/* + * 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 raft::linalg::detail { + +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); + + detail::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) +{ + detail::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 raft::linalg::detail diff --git a/cpp/include/raft/linalg/pca.cuh b/cpp/include/raft/linalg/pca.cuh new file mode 100644 index 0000000000..1384b13d3d --- /dev/null +++ b/cpp/include/raft/linalg/pca.cuh @@ -0,0 +1,194 @@ +/* + * SPDX-FileCopyrightText: Copyright (c) 2018-2026, NVIDIA CORPORATION. + * SPDX-License-Identifier: Apache-2.0 + */ + +#pragma once + +#include "detail/pca.cuh" + +namespace raft::linalg { + +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) +{ + detail::truncCompExpVars( + handle, in, components, explained_var, explained_var_ratio, noise_vars, prms, stream); +} + +/** + * @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) +{ + detail::pcaFit(handle, + input, + components, + explained_var, + explained_var_ratio, + singular_vals, + mu, + noise_vars, + prms, + stream, + 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) +{ + detail::pcaFitTransform(handle, + input, + trans_input, + components, + explained_var, + explained_var_ratio, + singular_vals, + mu, + noise_vars, + prms, + stream, + flip_signs_based_on_U); +} + +// TODO: implement pcaGetCovariance function +template +void pcaGetCovariance() +{ + detail::pcaGetCovariance(); +} + +// TODO: implement pcaGetPrecision function +template +void pcaGetPrecision() +{ + detail::pcaGetPrecision(); +} + +/** + * @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) +{ + detail::pcaInverseTransform( + handle, trans_input, components, singular_vals, mu, input, prms, stream); +} + +// TODO: implement pcaScore function +template +void pcaScore() +{ + detail::pcaScore(); +} + +// TODO: implement pcaScoreSamples function +template +void pcaScoreSamples() +{ + detail::pcaScoreSamples(); +} + +/** + * @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) +{ + detail::pcaTransform(handle, input, components, trans_input, singular_vals, mu, prms, stream); +} + +}; // end namespace raft::linalg diff --git a/cpp/include/raft/linalg/pca_types.hpp b/cpp/include/raft/linalg/pca_types.hpp new file mode 100644 index 0000000000..ec70b1a267 --- /dev/null +++ b/cpp/include/raft/linalg/pca_types.hpp @@ -0,0 +1,78 @@ +/* + * SPDX-FileCopyrightText: Copyright (c) 2018-2026, NVIDIA CORPORATION. + * SPDX-License-Identifier: Apache-2.0 + */ + +#pragma once + +#include + +namespace raft::linalg { + +/** + * @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, +}; + +class params { + public: + std::size_t n_rows; + std::size_t n_cols; + int gpu_id = 0; +}; + +class paramsSolver : public params { + public: + // math_t tol = 0.0; + float tol = 0.0; + std::uint32_t n_iterations = 15; + int verbose = 0; +}; + +template +class paramsTSVDTemplate : public paramsSolver { + public: + std::size_t n_components = 1; + enum_solver algorithm = enum_solver::COV_EIG_DQ; +}; + +/** + * @brief structure for pca parameters. Ref: + * http://scikit-learn.org/stable/modules/generated/sklearn.decomposition.PCA.html + * @param n_components: Number of components to keep. if n_components is not set all components are + * kept: + * @param copy: 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. + * @param whiten: 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. + * @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 verbose: 0: no error message printing, 1: print error messages + */ + +template +class paramsPCATemplate : public paramsTSVDTemplate { + public: + bool copy = true; // TODO unused, see #2830 and #2833 + bool whiten = false; +}; + +typedef paramsTSVDTemplate<> paramsTSVD; +typedef paramsPCATemplate<> paramsPCA; + +enum class mg_solver { COV_EIG_DQ, COV_EIG_JACOBI }; + +typedef paramsPCATemplate paramsPCAMG; +typedef paramsTSVDTemplate paramsTSVDMG; + +}; // end namespace raft::linalg diff --git a/cpp/include/raft/linalg/tsvd.cuh b/cpp/include/raft/linalg/tsvd.cuh new file mode 100644 index 0000000000..bb01c6e65d --- /dev/null +++ b/cpp/include/raft/linalg/tsvd.cuh @@ -0,0 +1,199 @@ +/* + * SPDX-FileCopyrightText: Copyright (c) 2018-2026, NVIDIA CORPORATION. + * SPDX-License-Identifier: Apache-2.0 + */ + +#pragma once + +#include "detail/tsvd.cuh" + +namespace raft::linalg { + +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) +{ + detail::calCompExpVarsSvd( + handle, in, components, singular_vals, explained_vars, explained_var_ratio, prms, 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) +{ + detail::calEig(handle, in, components, explained_var, prms, stream); +} + +/** + * @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) +{ + detail::signFlipComponents(handle, + input, + components, + n_samples, + n_features, + n_components, + stream, + center, + flip_signs_based_on_U); +} + +/** + * @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) +{ + detail::signFlip(input, n_rows, n_cols, components, n_cols_comp, stream); +} + +/** + * @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) +{ + detail::tsvdFit(handle, input, components, singular_vals, prms, stream, 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) +{ + detail::tsvdFitTransform(handle, + input, + trans_input, + components, + explained_var, + explained_var_ratio, + singular_vals, + prms, + stream, + flip_signs_based_on_U); +} + +/** + * @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) +{ + detail::tsvdTransform(handle, input, components, trans_input, prms, 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) +{ + detail::tsvdInverseTransform(handle, trans_input, components, input, prms, stream); +} + +}; // end namespace raft::linalg diff --git a/cpp/tests/CMakeLists.txt b/cpp/tests/CMakeLists.txt index dff227dd8b..b5c20e0212 100644 --- a/cpp/tests/CMakeLists.txt +++ b/cpp/tests/CMakeLists.txt @@ -1,6 +1,6 @@ # ============================================================================= # cmake-format: off -# SPDX-FileCopyrightText: Copyright (c) 2021-2025, NVIDIA CORPORATION. +# SPDX-FileCopyrightText: Copyright (c) 2021-2026, NVIDIA CORPORATION. # SPDX-License-Identifier: Apache-2.0 # cmake-format: on # ============================================================================= @@ -162,6 +162,7 @@ if(BUILD_TESTS) linalg/multiply.cu linalg/norm.cu linalg/normalize.cu + linalg/pca.cu linalg/power.cu linalg/randomized_svd.cu linalg/reduce.cu @@ -174,6 +175,7 @@ if(BUILD_TESTS) linalg/svd.cu linalg/ternary_op.cu linalg/transpose.cu + linalg/tsvd.cu linalg/unary_op.cu GPUS 1 diff --git a/cpp/tests/linalg/pca.cu b/cpp/tests/linalg/pca.cu new file mode 100644 index 0000000000..a3256a65aa --- /dev/null +++ b/cpp/tests/linalg/pca.cu @@ -0,0 +1,316 @@ +/* + * 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 + +namespace raft::linalg { + +template +struct PcaInputs { + T tolerance; + int len; + int n_row; + int n_col; + int len2; + 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; +} + +template +class PcaTest : public ::testing::TestWithParam> { + public: + PcaTest() + : params(::testing::TestWithParam>::GetParam()), + stream(handle.get_stream()), + 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.len, stream), + trans_data_ref(params.len, stream), + data(params.len, stream), + data_back(params.len, stream), + data2(params.len2, stream), + data2_back(params.len2, stream) + { + basicTest(); + advancedTest(); + } + + protected: + void basicTest() + { + raft::random::Rng r(params.seed, raft::random::GenPC); + int len = params.len; + + 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; + rmm::device_uvector explained_var_ratio(params.n_col, stream); + rmm::device_uvector singular_vals(params.n_col, stream); + rmm::device_uvector mean(params.n_col, stream); + rmm::device_uvector noise_vars(1, stream); + + 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); + + paramsPCA prms; + prms.n_cols = params.n_col; + prms.n_rows = params.n_row; + prms.n_components = params.n_col; + prms.whiten = false; + if (params.algo == 0) + prms.algorithm = solver::COV_EIG_DQ; + else + prms.algorithm = solver::COV_EIG_JACOBI; + + pcaFit(handle, + data.data(), + components.data(), + explained_vars.data(), + explained_var_ratio.data(), + singular_vals.data(), + mean.data(), + noise_vars.data(), + prms, + stream); + pcaTransform(handle, + data.data(), + components.data(), + trans_data.data(), + singular_vals.data(), + mean.data(), + prms, + stream); + pcaInverseTransform(handle, + trans_data.data(), + components.data(), + singular_vals.data(), + mean.data(), + data_back.data(), + prms, + stream); + } + + void advancedTest() + { + raft::random::Rng r(params.seed, raft::random::GenPC); + int len = params.len2; + + paramsPCA prms; + prms.n_cols = params.n_col2; + prms.n_rows = params.n_row2; + prms.n_components = params.n_col2; + prms.whiten = false; + if (params.algo == 0) + prms.algorithm = solver::COV_EIG_DQ; + else if (params.algo == 1) + prms.algorithm = solver::COV_EIG_JACOBI; + + r.uniform(data2.data(), len, T(-1.0), T(1.0), stream); + rmm::device_uvector data2_trans(prms.n_rows * prms.n_components, stream); + + int len_comp = params.n_col2 * prms.n_components; + rmm::device_uvector components2(len_comp, stream); + rmm::device_uvector explained_vars2(prms.n_components, stream); + rmm::device_uvector explained_var_ratio2(prms.n_components, stream); + rmm::device_uvector singular_vals2(prms.n_components, stream); + rmm::device_uvector mean2(prms.n_cols, stream); + rmm::device_uvector noise_vars2(1, stream); + + pcaFitTransform(handle, + data2.data(), + data2_trans.data(), + components2.data(), + explained_vars2.data(), + explained_var_ratio2.data(), + singular_vals2.data(), + mean2.data(), + noise_vars2.data(), + prms, + stream); + + pcaInverseTransform(handle, + data2_trans.data(), + components2.data(), + singular_vals2.data(), + mean2.data(), + data2_back.data(), + prms, + stream); + } + + protected: + raft::handle_t handle; + cudaStream_t stream = 0; + + PcaInputs params; + + 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, 3, 2, 1024 * 128, 1024, 128, 1234ULL, 0}, + {0.01f, 3 * 2, 3, 2, 256 * 32, 256, 32, 1234ULL, 1}}; + +const std::vector> inputsd2 = { + {0.01, 3 * 2, 3, 2, 1024 * 128, 1024, 128, 1234ULL, 0}, + {0.01, 3 * 2, 3, 2, 256 * 32, 256, 32, 1234ULL, 1}}; + +typedef PcaTest PcaTestValF; +TEST_P(PcaTestValF, Result) +{ + ASSERT_TRUE(devArrMatch(explained_vars.data(), + explained_vars_ref.data(), + params.n_col, + raft::CompareApprox(params.tolerance), + handle.get_stream())); +} + +typedef PcaTest PcaTestValD; +TEST_P(PcaTestValD, Result) +{ + ASSERT_TRUE(devArrMatch(explained_vars.data(), + explained_vars_ref.data(), + params.n_col, + raft::CompareApprox(params.tolerance), + handle.get_stream())); +} + +typedef PcaTest PcaTestLeftVecF; +TEST_P(PcaTestLeftVecF, Result) +{ + ASSERT_TRUE(devArrMatch(components.data(), + components_ref.data(), + (params.n_col * params.n_col), + raft::CompareApprox(params.tolerance), + handle.get_stream())); +} + +typedef PcaTest PcaTestLeftVecD; +TEST_P(PcaTestLeftVecD, Result) +{ + ASSERT_TRUE(devArrMatch(components.data(), + components_ref.data(), + (params.n_col * params.n_col), + raft::CompareApprox(params.tolerance), + handle.get_stream())); +} + +typedef PcaTest PcaTestTransDataF; +TEST_P(PcaTestTransDataF, Result) +{ + ASSERT_TRUE(devArrMatch(trans_data.data(), + trans_data_ref.data(), + (params.n_row * params.n_col), + raft::CompareApprox(params.tolerance), + handle.get_stream())); +} + +typedef PcaTest PcaTestTransDataD; +TEST_P(PcaTestTransDataD, Result) +{ + ASSERT_TRUE(devArrMatch(trans_data.data(), + trans_data_ref.data(), + (params.n_row * params.n_col), + raft::CompareApprox(params.tolerance), + handle.get_stream())); +} + +typedef PcaTest PcaTestDataVecSmallF; +TEST_P(PcaTestDataVecSmallF, Result) +{ + ASSERT_TRUE(devArrMatch(data.data(), + data_back.data(), + (params.n_col * params.n_col), + raft::CompareApprox(params.tolerance), + handle.get_stream())); +} + +typedef PcaTest PcaTestDataVecSmallD; +TEST_P(PcaTestDataVecSmallD, Result) +{ + ASSERT_TRUE(devArrMatch(data.data(), + data_back.data(), + (params.n_col * params.n_col), + raft::CompareApprox(params.tolerance), + handle.get_stream())); +} + +// FIXME: These tests are disabled due to driver 418+ making them fail: +// https://github.com/rapidsai/cuml/issues/379 +typedef PcaTest PcaTestDataVecF; +TEST_P(PcaTestDataVecF, Result) +{ + ASSERT_TRUE(devArrMatch(data2.data(), + data2_back.data(), + (params.n_col2 * params.n_col2), + raft::CompareApprox(params.tolerance), + handle.get_stream())); +} + +typedef PcaTest PcaTestDataVecD; +TEST_P(PcaTestDataVecD, Result) +{ + ASSERT_TRUE(devArrMatch(data2.data(), + data2_back.data(), + (params.n_col2 * params.n_col2), + raft::CompareApprox(params.tolerance), + handle.get_stream())); +} + +INSTANTIATE_TEST_CASE_P(PcaTests, PcaTestValF, ::testing::ValuesIn(inputsf2)); + +INSTANTIATE_TEST_CASE_P(PcaTests, PcaTestValD, ::testing::ValuesIn(inputsd2)); + +INSTANTIATE_TEST_CASE_P(PcaTests, PcaTestLeftVecF, ::testing::ValuesIn(inputsf2)); + +INSTANTIATE_TEST_CASE_P(PcaTests, PcaTestLeftVecD, ::testing::ValuesIn(inputsd2)); + +INSTANTIATE_TEST_CASE_P(PcaTests, PcaTestDataVecSmallF, ::testing::ValuesIn(inputsf2)); + +INSTANTIATE_TEST_CASE_P(PcaTests, PcaTestDataVecSmallD, ::testing::ValuesIn(inputsd2)); + +INSTANTIATE_TEST_CASE_P(PcaTests, PcaTestTransDataF, ::testing::ValuesIn(inputsf2)); + +INSTANTIATE_TEST_CASE_P(PcaTests, PcaTestTransDataD, ::testing::ValuesIn(inputsd2)); + +INSTANTIATE_TEST_CASE_P(PcaTests, PcaTestDataVecF, ::testing::ValuesIn(inputsf2)); + +INSTANTIATE_TEST_CASE_P(PcaTests, PcaTestDataVecD, ::testing::ValuesIn(inputsd2)); + +} // end namespace raft::linalg diff --git a/cpp/tests/linalg/tsvd.cu b/cpp/tests/linalg/tsvd.cu new file mode 100644 index 0000000000..6892fade3c --- /dev/null +++ b/cpp/tests/linalg/tsvd.cu @@ -0,0 +1,206 @@ +/* + * 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 + +namespace raft::linalg { + +template +struct TsvdInputs { + T tolerance; + int n_row; + int n_col; + int n_row2; + int n_col2; + float redundancy; + unsigned long long int seed; + int algo; +}; + +template +::std::ostream& operator<<(::std::ostream& os, const TsvdInputs& dims) +{ + return os; +} + +template +class TsvdTest : public ::testing::TestWithParam> { + public: + TsvdTest() + : params(::testing::TestWithParam>::GetParam()), + stream(handle.get_stream()), + components(0, stream), + components_ref(0, stream), + data2(0, stream), + data2_back(0, stream) + { + basicTest(); + advancedTest(); + } + + protected: + void basicTest() + { + raft::random::Rng r(params.seed, raft::random::GenPC); + int len = params.n_row * params.n_col; + + rmm::device_uvector data(len, stream); + + std::vector data_h = {1.0, 2.0, 4.0, 2.0, 4.0, 5.0, 5.0, 4.0, 2.0, 1.0, 6.0, 4.0}; + data_h.resize(len); + raft::update_device(data.data(), data_h.data(), len, stream); + + int len_comp = params.n_col * params.n_col; + components.resize(len_comp, stream); + rmm::device_uvector singular_vals(params.n_col, stream); + + std::vector components_ref_h = { + 0.3951, 0.1532, 0.9058, 0.7111, -0.6752, -0.1959, 0.5816, 0.7215, -0.3757}; + components_ref_h.resize(len_comp); + + components_ref.resize(len_comp, stream); + raft::update_device(components_ref.data(), components_ref_h.data(), len_comp, stream); + + paramsTSVD prms; + prms.n_cols = params.n_col; + prms.n_rows = params.n_row; + prms.n_components = params.n_col; + if (params.algo == 0) + prms.algorithm = solver::COV_EIG_DQ; + else + prms.algorithm = solver::COV_EIG_JACOBI; + + tsvdFit(handle, data.data(), components.data(), singular_vals.data(), prms, stream); + } + + void advancedTest() + { + raft::random::Rng r(params.seed, raft::random::GenPC); + int len = params.n_row2 * params.n_col2; + + paramsTSVD prms; + prms.n_cols = params.n_col2; + prms.n_rows = params.n_row2; + prms.n_components = params.n_col2; + if (params.algo == 0) + prms.algorithm = solver::COV_EIG_DQ; + else if (params.algo == 1) + prms.algorithm = solver::COV_EIG_JACOBI; + else + prms.n_components = params.n_col2 - 15; + + data2.resize(len, stream); + int redundant_cols = int(params.redundancy * params.n_col2); + int redundant_len = params.n_row2 * redundant_cols; + + int informative_cols = params.n_col2 - redundant_cols; + int informative_len = params.n_row2 * informative_cols; + + r.uniform(data2.data(), informative_len, T(-1.0), T(1.0), stream); + RAFT_CUDA_TRY(cudaMemcpyAsync(data2.data() + informative_len, + data2.data(), + redundant_len * sizeof(T), + cudaMemcpyDeviceToDevice, + stream)); + rmm::device_uvector data2_trans(prms.n_rows * prms.n_components, stream); + + int len_comp = params.n_col2 * prms.n_components; + rmm::device_uvector components2(len_comp, stream); + rmm::device_uvector explained_vars2(prms.n_components, stream); + rmm::device_uvector explained_var_ratio2(prms.n_components, stream); + rmm::device_uvector singular_vals2(prms.n_components, stream); + + tsvdFitTransform(handle, + data2.data(), + data2_trans.data(), + components2.data(), + explained_vars2.data(), + explained_var_ratio2.data(), + singular_vals2.data(), + prms, + stream); + + data2_back.resize(len, stream); + tsvdInverseTransform( + handle, data2_trans.data(), components2.data(), data2_back.data(), prms, stream); + } + + protected: + raft::handle_t handle; + cudaStream_t stream = 0; + + TsvdInputs params; + rmm::device_uvector components, components_ref, data2, data2_back; +}; + +const std::vector> inputsf2 = {{0.01f, 4, 3, 1024, 128, 0.25f, 1234ULL, 0}, + {0.01f, 4, 3, 1024, 128, 0.25f, 1234ULL, 1}, + {0.04f, 4, 3, 512, 64, 0.25f, 1234ULL, 2}, + {0.04f, 4, 3, 512, 64, 0.25f, 1234ULL, 2}}; + +const std::vector> inputsd2 = {{0.01, 4, 3, 1024, 128, 0.25f, 1234ULL, 0}, + {0.01, 4, 3, 1024, 128, 0.25f, 1234ULL, 1}, + {0.05, 4, 3, 512, 64, 0.25f, 1234ULL, 2}, + {0.05, 4, 3, 512, 64, 0.25f, 1234ULL, 2}}; + +typedef TsvdTest TsvdTestLeftVecF; +TEST_P(TsvdTestLeftVecF, Result) +{ + ASSERT_TRUE(devArrMatch(components.data(), + components_ref.data(), + (params.n_col * params.n_col), + raft::CompareApprox(params.tolerance), + handle.get_stream())); +} + +typedef TsvdTest TsvdTestLeftVecD; +TEST_P(TsvdTestLeftVecD, Result) +{ + ASSERT_TRUE(devArrMatch(components.data(), + components_ref.data(), + (params.n_col * params.n_col), + raft::CompareApprox(params.tolerance), + handle.get_stream())); +} + +typedef TsvdTest TsvdTestDataVecF; +TEST_P(TsvdTestDataVecF, Result) +{ + ASSERT_TRUE(devArrMatch(data2.data(), + data2_back.data(), + (params.n_col2 * params.n_col2), + raft::CompareApprox(params.tolerance), + handle.get_stream())); +} + +typedef TsvdTest TsvdTestDataVecD; +TEST_P(TsvdTestDataVecD, Result) +{ + ASSERT_TRUE(devArrMatch(data2.data(), + data2_back.data(), + (params.n_col2 * params.n_col2), + raft::CompareApprox(params.tolerance), + handle.get_stream())); +} + +INSTANTIATE_TEST_CASE_P(TsvdTests, TsvdTestLeftVecF, ::testing::ValuesIn(inputsf2)); + +INSTANTIATE_TEST_CASE_P(TsvdTests, TsvdTestLeftVecD, ::testing::ValuesIn(inputsd2)); + +INSTANTIATE_TEST_CASE_P(TsvdTests, TsvdTestDataVecF, ::testing::ValuesIn(inputsf2)); + +INSTANTIATE_TEST_CASE_P(TsvdTests, TsvdTestDataVecD, ::testing::ValuesIn(inputsd2)); + +} // end namespace raft::linalg From 7289840b412e4fa4502d9f2ab182775742504354 Mon Sep 17 00:00:00 2001 From: aamijar Date: Sat, 14 Feb 2026 01:50:11 +0000 Subject: [PATCH 02/12] mdspan public api --- cpp/include/raft/linalg/detail/pca.cuh | 192 +++++++-------- cpp/include/raft/linalg/detail/tsvd.cuh | 132 +++++----- cpp/include/raft/linalg/pca.cuh | 312 ++++++++++++------------ cpp/include/raft/linalg/tsvd.cuh | 285 ++++++++++------------ cpp/tests/linalg/pca.cu | 139 ++++++----- cpp/tests/linalg/tsvd.cu | 60 +++-- 6 files changed, 556 insertions(+), 564 deletions(-) diff --git a/cpp/include/raft/linalg/detail/pca.cuh b/cpp/include/raft/linalg/detail/pca.cuh index 43d7640efd..eff409c730 100644 --- a/cpp/include/raft/linalg/detail/pca.cuh +++ b/cpp/include/raft/linalg/detail/pca.cuh @@ -5,7 +5,7 @@ #pragma once -#include +#include #include #include #include @@ -26,7 +26,7 @@ namespace raft::linalg::detail { template -void truncCompExpVars(const raft::handle_t& handle, +void truncCompExpVars(raft::resources const& handle, math_t* in, math_t* components, math_t* explained_var, @@ -101,7 +101,7 @@ void truncCompExpVars(const raft::handle_t& handle, * @param[in] stream cuda stream */ template -void pcaFit(const raft::handle_t& handle, +void pcaFit(raft::resources const& handle, math_t* input, math_t* components, math_t* explained_var, @@ -113,7 +113,7 @@ void pcaFit(const raft::handle_t& handle, cudaStream_t stream, bool flip_signs_based_on_U = false) { - auto cublas_handle = handle.get_cublas_handle(); + auto cublas_handle = raft::resource::get_cublas_handle(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"); @@ -157,63 +157,52 @@ void pcaFit(const raft::handle_t& handle, } /** - * @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). + * @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 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) +void pcaTransform(raft::resources const& handle, + math_t* input, + math_t* components, + math_t* trans_input, + math_t* singular_vals, + math_t* mu, + const paramsPCA& prms, + cudaStream_t stream) { - detail::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); -} + 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"); -// TODO: implement pcaGetCovariance function -template -void pcaGetCovariance() -{ - ASSERT(false, "pcaGetCovariance: will be implemented!"); -} + 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); -// TODO: implement pcaGetPrecision function -template -void pcaGetPrecision() -{ - ASSERT(false, "pcaGetPrecision: will be implemented!"); + 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); + detail::tsvdTransform(handle, input, components_copy.data(), trans_input, prms, stream); + raft::stats::meanAdd(input, input, mu, prms.n_cols, prms.n_rows, stream); } /** @@ -230,7 +219,7 @@ void pcaGetPrecision() * @param[in] stream cuda stream */ template -void pcaInverseTransform(const raft::handle_t& handle, +void pcaInverseTransform(raft::resources const& handle, math_t* trans_input, math_t* components, math_t* singular_vals, @@ -267,67 +256,50 @@ void pcaInverseTransform(const raft::handle_t& handle, 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 + * @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 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) +void pcaFitTransform(raft::resources const& 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) { - 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); - detail::tsvdTransform(handle, input, components_copy.data(), trans_input, prms, stream); - raft::stats::meanAdd(input, input, mu, prms.n_cols, prms.n_rows, stream); + detail::pcaFit(handle, + input, + components, + explained_var, + explained_var_ratio, + singular_vals, + mu, + noise_vars, + prms, + stream, + flip_signs_based_on_U); + detail::pcaTransform(handle, input, components, trans_input, singular_vals, mu, prms, stream); } }; // end namespace raft::linalg::detail diff --git a/cpp/include/raft/linalg/detail/tsvd.cuh b/cpp/include/raft/linalg/detail/tsvd.cuh index 39adffc8e3..8dd95cee1c 100644 --- a/cpp/include/raft/linalg/detail/tsvd.cuh +++ b/cpp/include/raft/linalg/detail/tsvd.cuh @@ -5,7 +5,7 @@ #pragma once -#include +#include #include #include #include @@ -38,7 +38,7 @@ namespace raft::linalg::detail { template -void calCompExpVarsSvd(const raft::handle_t& handle, +void calCompExpVarsSvd(raft::resources const& handle, math_t* in, math_t* components, math_t* singular_vals, @@ -47,8 +47,8 @@ void calCompExpVarsSvd(const raft::handle_t& handle, const paramsTSVD& prms, cudaStream_t stream) { - auto cusolver_handle = handle.get_cusolver_dn_handle(); - auto cublas_handle = handle.get_cublas_handle(); + auto cusolver_handle = raft::resource::get_cusolver_dn_handle(handle); + auto cublas_handle = raft::resource::get_cublas_handle(handle); auto diff = prms.n_cols - prms.n_components; math_t ratio = math_t(diff) / math_t(prms.n_cols); @@ -97,14 +97,14 @@ void calCompExpVarsSvd(const raft::handle_t& handle, } template -void calEig(const raft::handle_t& handle, +void calEig(raft::resources const& 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(); + auto cusolver_handle = raft::resource::get_cusolver_dn_handle(handle); if (prms.algorithm == enum_solver::COV_EIG_JACOBI) { raft::linalg::eigJacobi(handle, @@ -146,7 +146,7 @@ void calEig(const raft::handle_t& handle, * @{ */ template -void signFlipComponents(const raft::handle_t& handle, +void signFlipComponents(raft::resources const& handle, math_t* input, math_t* components, std::size_t n_samples, @@ -299,7 +299,7 @@ void signFlip(math_t* input, * @param[in] stream cuda stream */ template -void tsvdFit(const raft::handle_t& handle, +void tsvdFit(raft::resources const& handle, math_t* input, math_t* components, math_t* singular_vals, @@ -307,7 +307,7 @@ void tsvdFit(const raft::handle_t& handle, cudaStream_t stream, bool flip_signs_based_on_U = false) { - auto cublas_handle = handle.get_cublas_handle(); + auto cublas_handle = raft::resource::get_cublas_handle(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"); @@ -369,61 +369,6 @@ void tsvdFit(const raft::handle_t& handle, 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) -{ - detail::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 @@ -434,7 +379,7 @@ void tsvdFitTransform(const raft::handle_t& handle, * @param[in] stream cuda stream */ template -void tsvdTransform(const raft::handle_t& handle, +void tsvdTransform(raft::resources const& handle, math_t* input, math_t* components, math_t* trans_input, @@ -475,7 +420,7 @@ void tsvdTransform(const raft::handle_t& handle, * @param[in] stream cuda stream */ template -void tsvdInverseTransform(const raft::handle_t& handle, +void tsvdInverseTransform(raft::resources const& handle, math_t* trans_input, math_t* components, math_t* input, @@ -505,4 +450,59 @@ void tsvdInverseTransform(const raft::handle_t& handle, stream); } +/** + * @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(raft::resources const& 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) +{ + detail::tsvdFit(handle, input, components, singular_vals, prms, stream, flip_signs_based_on_U); + detail::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); + raft::resource::sync_stream(handle, stream); + math_t scalar = math_t(1) / total_vars_h; + + raft::linalg::scalarMultiply( + explained_var_ratio, explained_var, scalar, prms.n_components, stream); +} + }; // end namespace raft::linalg::detail diff --git a/cpp/include/raft/linalg/pca.cuh b/cpp/include/raft/linalg/pca.cuh index 1384b13d3d..9c305c75ff 100644 --- a/cpp/include/raft/linalg/pca.cuh +++ b/cpp/include/raft/linalg/pca.cuh @@ -7,188 +7,198 @@ #include "detail/pca.cuh" +#include +#include + namespace raft::linalg { -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) -{ - detail::truncCompExpVars( - handle, in, components, explained_var, explained_var_ratio, noise_vars, prms, stream); -} +/** + * @defgroup pca PCA operations + * @{ + */ /** - * @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 + * @brief perform fit operation for PCA. Generates eigenvectors, explained vars, singular vals, etc. + * @tparam math_t data-type upon which the math operation will be performed + * @tparam idx_t integer type used for indexing + * @param[in] handle raft::resources + * @param[in] prms PCA parameters (n_components, algorithm, whiten, etc.) + * @param[inout] input the data is fitted to PCA. Size n_rows x n_cols (col-major). Modified + * temporarily during computation. + * @param[out] components the principal components of the input data. Size n_components x n_cols + * (col-major). + * @param[out] explained_var explained variances (eigenvalues) of the principal components. Size + * n_components. + * @param[out] explained_var_ratio the ratio of the explained variance and total variance. Size + * n_components. + * @param[out] singular_vals singular values of the data. Size n_components. + * @param[out] mu mean of all the features (all the columns in the data). Size n_cols. + * @param[out] noise_vars variance of the noise. Scalar. + * @param[in] flip_signs_based_on_U whether to determine signs by U (true) or V.T (false) */ -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) +template +void pca_fit(raft::resources const& handle, + const paramsPCA& prms, + 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) { + auto stream = resource::get_cuda_stream(handle); + + paramsPCA prms_with_dims = prms; + prms_with_dims.n_rows = static_cast(input.extent(0)); + prms_with_dims.n_cols = static_cast(input.extent(1)); + detail::pcaFit(handle, - input, - components, - explained_var, - explained_var_ratio, - singular_vals, - mu, - noise_vars, - prms, + input.data_handle(), + components.data_handle(), + explained_var.data_handle(), + explained_var_ratio.data_handle(), + singular_vals.data_handle(), + mu.data_handle(), + noise_vars.data_handle(), + prms_with_dims, stream, flip_signs_based_on_U); } /** - * @brief perform fit and transform operations for the pca. Generates transformed data, + * @brief perform fit and transform operations for 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 + * @tparam math_t data-type upon which the math operation will be performed + * @tparam idx_t integer type used for indexing + * @param[in] handle raft::resources + * @param[in] prms PCA parameters (n_components, algorithm, whiten, etc.) + * @param[inout] input the data is fitted to PCA. Size n_rows x n_cols (col-major). Modified + * temporarily during computation. + * @param[out] trans_input the transformed data. Size n_rows x n_components (col-major). + * @param[out] components the principal components of the input data. Size n_components x n_cols + * (col-major). + * @param[out] explained_var explained variances (eigenvalues) of the principal components. Size + * n_components. + * @param[out] explained_var_ratio the ratio of the explained variance and total variance. Size + * n_components. + * @param[out] singular_vals singular values of the data. Size n_components. + * @param[out] mu mean of all the features (all the columns in the data). Size n_cols. + * @param[out] noise_vars variance of the noise. Scalar. + * @param[in] flip_signs_based_on_U whether to determine signs by U (true) or V.T (false) */ -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) +template +void pca_fit_transform(raft::resources const& handle, + const paramsPCA& prms, + 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) { + auto stream = resource::get_cuda_stream(handle); + + paramsPCA prms_with_dims = prms; + prms_with_dims.n_rows = static_cast(input.extent(0)); + prms_with_dims.n_cols = static_cast(input.extent(1)); + detail::pcaFitTransform(handle, - input, - trans_input, - components, - explained_var, - explained_var_ratio, - singular_vals, - mu, - noise_vars, - prms, + input.data_handle(), + trans_input.data_handle(), + components.data_handle(), + explained_var.data_handle(), + explained_var_ratio.data_handle(), + singular_vals.data_handle(), + mu.data_handle(), + noise_vars.data_handle(), + prms_with_dims, stream, flip_signs_based_on_U); } -// TODO: implement pcaGetCovariance function -template -void pcaGetCovariance() -{ - detail::pcaGetCovariance(); -} - -// TODO: implement pcaGetPrecision function -template -void pcaGetPrecision() -{ - detail::pcaGetPrecision(); -} - /** - * @brief performs inverse transform operation for the pca. Transforms the transformed data back to + * @brief performs inverse transform operation for 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 + * @tparam math_t data-type upon which the math operation will be performed + * @tparam idx_t integer type used for indexing + * @param[in] handle raft::resources + * @param[in] prms PCA parameters (n_components, algorithm, whiten, etc.) + * @param[in] trans_input the transformed data. Size n_rows x n_components (col-major). + * @param[in] components the principal components of the input data. Size n_components x n_cols + * (col-major). + * @param[in] singular_vals singular values of the data. Size n_components. + * @param[in] mu mean of features (every column). Size n_cols. + * @param[out] output the reconstructed data. Size n_rows x n_cols (col-major). */ -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) +template +void pca_inverse_transform(raft::resources const& handle, + const paramsPCA& prms, + 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::pcaInverseTransform( - handle, trans_input, components, singular_vals, mu, input, prms, stream); -} + auto stream = resource::get_cuda_stream(handle); -// TODO: implement pcaScore function -template -void pcaScore() -{ - detail::pcaScore(); -} + paramsPCA prms_with_dims = prms; + prms_with_dims.n_rows = static_cast(output.extent(0)); + prms_with_dims.n_cols = static_cast(output.extent(1)); -// TODO: implement pcaScoreSamples function -template -void pcaScoreSamples() -{ - detail::pcaScoreSamples(); + detail::pcaInverseTransform(handle, + trans_input.data_handle(), + components.data_handle(), + singular_vals.data_handle(), + mu.data_handle(), + output.data_handle(), + prms_with_dims, + stream); } /** - * @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 + * @brief performs transform operation for PCA. Transforms the data to eigenspace. + * @tparam math_t data-type upon which the math operation will be performed + * @tparam idx_t integer type used for indexing + * @param[in] handle raft::resources + * @param[in] prms PCA parameters (n_components, algorithm, whiten, etc.) + * @param[inout] input the data to be transformed. Size n_rows x n_cols (col-major). Modified + * temporarily during computation (mean-centered then restored). + * @param[in] components principal components of the input data. Size n_components x n_cols + * (col-major). + * @param[in] singular_vals singular values of the data. Size n_components. + * @param[in] mu mean value of the input data. Size n_cols. + * @param[out] trans_input the transformed data. Size n_rows x n_components (col-major). */ -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) +template +void pca_transform(raft::resources const& handle, + const paramsPCA& prms, + 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::pcaTransform(handle, input, components, trans_input, singular_vals, mu, prms, stream); + auto stream = resource::get_cuda_stream(handle); + + paramsPCA prms_with_dims = prms; + prms_with_dims.n_rows = static_cast(input.extent(0)); + prms_with_dims.n_cols = static_cast(input.extent(1)); + + detail::pcaTransform(handle, + input.data_handle(), + components.data_handle(), + trans_input.data_handle(), + singular_vals.data_handle(), + mu.data_handle(), + prms_with_dims, + stream); } +/** @} */ // end group pca + }; // end namespace raft::linalg diff --git a/cpp/include/raft/linalg/tsvd.cuh b/cpp/include/raft/linalg/tsvd.cuh index bb01c6e65d..8a54e4cf1b 100644 --- a/cpp/include/raft/linalg/tsvd.cuh +++ b/cpp/include/raft/linalg/tsvd.cuh @@ -7,193 +7,164 @@ #include "detail/tsvd.cuh" -namespace raft::linalg { - -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) -{ - detail::calCompExpVarsSvd( - handle, in, components, singular_vals, explained_vars, explained_var_ratio, prms, stream); -} +#include +#include -template -void calEig(const raft::handle_t& handle, - math_t* in, - math_t* components, - math_t* explained_var, - const paramsTSVDTemplate& prms, - cudaStream_t stream) -{ - detail::calEig(handle, in, components, explained_var, prms, stream); -} +namespace raft::linalg { /** - * @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) + * @defgroup tsvd Truncated SVD operations * @{ */ -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) -{ - detail::signFlipComponents(handle, - input, - components, - n_samples, - n_features, - n_components, - stream, - center, - flip_signs_based_on_U); -} /** - * @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 - * @{ + * @brief perform fit operation for tSVD. Generates eigenvectors, singular vals, etc. + * @tparam math_t data-type upon which the math operation will be performed + * @tparam idx_t integer type used for indexing + * @param[in] handle raft::resources + * @param[in] prms data structure that includes all the parameters from input size to algorithm. + * @param[inout] input the data is fitted to tSVD. Size n_rows x n_cols (col-major). + * @param[out] components the principal components of the input data. Size n_components x n_cols + * (col-major). + * @param[out] singular_vals singular values of the data. Size n_components. + * @param[in] flip_signs_based_on_U whether to determine signs by U (true) or V.T (false) */ -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) +template +void tsvd_fit(raft::resources const& handle, + const paramsTSVD& prms, + raft::device_matrix_view input, + raft::device_matrix_view components, + raft::device_vector_view singular_vals, + bool flip_signs_based_on_U = false) { - detail::signFlip(input, n_rows, n_cols, components, n_cols_comp, stream); -} + auto stream = resource::get_cuda_stream(handle); -/** - * @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) -{ - detail::tsvdFit(handle, input, components, singular_vals, prms, stream, flip_signs_based_on_U); + paramsTSVD prms_with_dims = prms; + prms_with_dims.n_rows = static_cast(input.extent(0)); + prms_with_dims.n_cols = static_cast(input.extent(1)); + + detail::tsvdFit(handle, + input.data_handle(), + components.data_handle(), + singular_vals.data_handle(), + prms_with_dims, + stream, + flip_signs_based_on_U); } /** - * @brief performs fit and transform operations for the tsvd. Generates transformed data, + * @brief performs fit and transform operations for 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 + * @tparam math_t data-type upon which the math operation will be performed + * @tparam idx_t integer type used for indexing + * @param[in] handle raft::resources + * @param[in] prms data structure that includes all the parameters from input size to algorithm. + * @param[inout] input the data is fitted to tSVD. Size n_rows x n_cols (col-major). + * @param[out] trans_input the transformed data. Size n_rows x n_components (col-major). + * @param[out] components the principal components of the input data. Size n_components x n_cols + * (col-major). + * @param[out] explained_var explained variances (eigenvalues) of the principal components. Size + * n_components. + * @param[out] explained_var_ratio the ratio of the explained variance and total variance. Size + * n_components. + * @param[out] singular_vals singular values of the data. Size n_components. + * @param[in] flip_signs_based_on_U whether to determine signs by U (true) or V.T (false) */ -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) +template +void tsvd_fit_transform(raft::resources const& handle, + const paramsTSVD& prms, + 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, + bool flip_signs_based_on_U = false) { + auto stream = resource::get_cuda_stream(handle); + + paramsTSVD prms_with_dims = prms; + prms_with_dims.n_rows = static_cast(input.extent(0)); + prms_with_dims.n_cols = static_cast(input.extent(1)); + detail::tsvdFitTransform(handle, - input, - trans_input, - components, - explained_var, - explained_var_ratio, - singular_vals, - prms, + input.data_handle(), + trans_input.data_handle(), + components.data_handle(), + explained_var.data_handle(), + explained_var_ratio.data_handle(), + singular_vals.data_handle(), + prms_with_dims, stream, flip_signs_based_on_U); } /** - * @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 + * @brief performs transform operation for tSVD. Transforms the data to eigenspace. + * @tparam math_t data-type upon which the math operation will be performed + * @tparam idx_t integer type used for indexing + * @param[in] handle raft::resources + * @param[in] prms data structure that includes all the parameters from input size to algorithm. + * @param[in] input the data to be transformed. Size n_rows x n_cols (col-major). + * @param[in] components principal components of the input data. Size n_components x n_cols + * (col-major). + * @param[out] trans_input output that is transformed version of input. Size n_rows x n_components + * (col-major). */ -template -void tsvdTransform(const raft::handle_t& handle, - math_t* input, - math_t* components, - math_t* trans_input, - const paramsTSVD& prms, - cudaStream_t stream) +template +void tsvd_transform(raft::resources const& handle, + const paramsTSVD& prms, + raft::device_matrix_view input, + raft::device_matrix_view components, + raft::device_matrix_view trans_input) { - detail::tsvdTransform(handle, input, components, trans_input, prms, stream); + auto stream = resource::get_cuda_stream(handle); + + paramsTSVD prms_with_dims = prms; + prms_with_dims.n_rows = static_cast(input.extent(0)); + prms_with_dims.n_cols = static_cast(input.extent(1)); + + detail::tsvdTransform(handle, + input.data_handle(), + components.data_handle(), + trans_input.data_handle(), + prms_with_dims, + stream); } /** - * @brief performs inverse transform operation for the tsvd. Transforms the transformed data back to + * @brief performs inverse transform operation for 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 + * @tparam math_t data-type upon which the math operation will be performed + * @tparam idx_t integer type used for indexing + * @param[in] handle raft::resources + * @param[in] prms data structure that includes all the parameters from input size to algorithm. + * @param[in] trans_input the transformed data. Size n_rows x n_components (col-major). + * @param[in] components transpose of the principal components. Size n_components x n_cols + * (col-major). + * @param[out] output the reconstructed data. Size n_rows x n_cols (col-major). */ -template -void tsvdInverseTransform(const raft::handle_t& handle, - math_t* trans_input, - math_t* components, - math_t* input, - const paramsTSVD& prms, - cudaStream_t stream) +template +void tsvd_inverse_transform(raft::resources const& handle, + const paramsTSVD& prms, + raft::device_matrix_view trans_input, + raft::device_matrix_view components, + raft::device_matrix_view output) { - detail::tsvdInverseTransform(handle, trans_input, components, input, prms, stream); + auto stream = resource::get_cuda_stream(handle); + + paramsTSVD prms_with_dims = prms; + prms_with_dims.n_rows = static_cast(output.extent(0)); + prms_with_dims.n_cols = static_cast(output.extent(1)); + + detail::tsvdInverseTransform(handle, + trans_input.data_handle(), + components.data_handle(), + output.data_handle(), + prms_with_dims, + stream); } +/** @} */ // end group tsvd + }; // end namespace raft::linalg diff --git a/cpp/tests/linalg/pca.cu b/cpp/tests/linalg/pca.cu index a3256a65aa..a0d87995fa 100644 --- a/cpp/tests/linalg/pca.cu +++ b/cpp/tests/linalg/pca.cu @@ -5,7 +5,8 @@ #include "../test_utils.cuh" -#include +#include +#include #include #include #include @@ -42,7 +43,7 @@ class PcaTest : public ::testing::TestWithParam> { public: PcaTest() : params(::testing::TestWithParam>::GetParam()), - stream(handle.get_stream()), + stream(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), @@ -97,32 +98,40 @@ class PcaTest : public ::testing::TestWithParam> { else prms.algorithm = solver::COV_EIG_JACOBI; - pcaFit(handle, - data.data(), - components.data(), - explained_vars.data(), - explained_var_ratio.data(), - singular_vals.data(), - mean.data(), - noise_vars.data(), - prms, - stream); - pcaTransform(handle, - data.data(), - components.data(), - trans_data.data(), - singular_vals.data(), - mean.data(), - prms, - stream); - pcaInverseTransform(handle, - trans_data.data(), - components.data(), - singular_vals.data(), - mean.data(), - data_back.data(), - prms, - stream); + auto input_view = raft::make_device_matrix_view( + data.data(), prms.n_rows, prms.n_cols); + auto components_view = raft::make_device_matrix_view( + components.data(), prms.n_components, prms.n_cols); + auto explained_var_view = + raft::make_device_vector_view(explained_vars.data(), prms.n_components); + auto explained_var_ratio_view = + raft::make_device_vector_view(explained_var_ratio.data(), prms.n_components); + auto singular_vals_view = + raft::make_device_vector_view(singular_vals.data(), prms.n_components); + auto mu_view = raft::make_device_vector_view(mean.data(), prms.n_cols); + auto noise_vars_view = raft::make_device_scalar_view(noise_vars.data()); + + pca_fit(handle, + prms, + input_view, + components_view, + explained_var_view, + explained_var_ratio_view, + singular_vals_view, + mu_view, + noise_vars_view); + + auto trans_data_view = raft::make_device_matrix_view( + trans_data.data(), prms.n_rows, prms.n_components); + + pca_transform( + handle, prms, input_view, components_view, singular_vals_view, mu_view, trans_data_view); + + auto data_back_view = raft::make_device_matrix_view( + data_back.data(), prms.n_rows, prms.n_cols); + + pca_inverse_transform( + handle, prms, trans_data_view, components_view, singular_vals_view, mu_view, data_back_view); } void advancedTest() @@ -151,31 +160,41 @@ class PcaTest : public ::testing::TestWithParam> { rmm::device_uvector mean2(prms.n_cols, stream); rmm::device_uvector noise_vars2(1, stream); - pcaFitTransform(handle, - data2.data(), - data2_trans.data(), - components2.data(), - explained_vars2.data(), - explained_var_ratio2.data(), - singular_vals2.data(), - mean2.data(), - noise_vars2.data(), - prms, - stream); - - pcaInverseTransform(handle, - data2_trans.data(), - components2.data(), - singular_vals2.data(), - mean2.data(), - data2_back.data(), - prms, - stream); + auto input_view = raft::make_device_matrix_view( + data2.data(), prms.n_rows, prms.n_cols); + auto trans_view = raft::make_device_matrix_view( + data2_trans.data(), prms.n_rows, prms.n_components); + auto comp_view = raft::make_device_matrix_view( + components2.data(), prms.n_components, prms.n_cols); + auto ev_view = + raft::make_device_vector_view(explained_vars2.data(), prms.n_components); + auto evr_view = + raft::make_device_vector_view(explained_var_ratio2.data(), prms.n_components); + auto sv_view = + raft::make_device_vector_view(singular_vals2.data(), prms.n_components); + auto mu_view = raft::make_device_vector_view(mean2.data(), prms.n_cols); + auto noise_view = raft::make_device_scalar_view(noise_vars2.data()); + + pca_fit_transform(handle, + prms, + input_view, + trans_view, + comp_view, + ev_view, + evr_view, + sv_view, + mu_view, + noise_view); + + auto data2_back_view = raft::make_device_matrix_view( + data2_back.data(), prms.n_rows, prms.n_cols); + + pca_inverse_transform(handle, prms, trans_view, comp_view, sv_view, mu_view, data2_back_view); } protected: - raft::handle_t handle; - cudaStream_t stream = 0; + raft::device_resources handle; + cudaStream_t stream; PcaInputs params; @@ -198,7 +217,7 @@ TEST_P(PcaTestValF, Result) explained_vars_ref.data(), params.n_col, raft::CompareApprox(params.tolerance), - handle.get_stream())); + resource::get_cuda_stream(handle))); } typedef PcaTest PcaTestValD; @@ -208,7 +227,7 @@ TEST_P(PcaTestValD, Result) explained_vars_ref.data(), params.n_col, raft::CompareApprox(params.tolerance), - handle.get_stream())); + resource::get_cuda_stream(handle))); } typedef PcaTest PcaTestLeftVecF; @@ -218,7 +237,7 @@ TEST_P(PcaTestLeftVecF, Result) components_ref.data(), (params.n_col * params.n_col), raft::CompareApprox(params.tolerance), - handle.get_stream())); + resource::get_cuda_stream(handle))); } typedef PcaTest PcaTestLeftVecD; @@ -228,7 +247,7 @@ TEST_P(PcaTestLeftVecD, Result) components_ref.data(), (params.n_col * params.n_col), raft::CompareApprox(params.tolerance), - handle.get_stream())); + resource::get_cuda_stream(handle))); } typedef PcaTest PcaTestTransDataF; @@ -238,7 +257,7 @@ TEST_P(PcaTestTransDataF, Result) trans_data_ref.data(), (params.n_row * params.n_col), raft::CompareApprox(params.tolerance), - handle.get_stream())); + resource::get_cuda_stream(handle))); } typedef PcaTest PcaTestTransDataD; @@ -248,7 +267,7 @@ TEST_P(PcaTestTransDataD, Result) trans_data_ref.data(), (params.n_row * params.n_col), raft::CompareApprox(params.tolerance), - handle.get_stream())); + resource::get_cuda_stream(handle))); } typedef PcaTest PcaTestDataVecSmallF; @@ -258,7 +277,7 @@ TEST_P(PcaTestDataVecSmallF, Result) data_back.data(), (params.n_col * params.n_col), raft::CompareApprox(params.tolerance), - handle.get_stream())); + resource::get_cuda_stream(handle))); } typedef PcaTest PcaTestDataVecSmallD; @@ -268,7 +287,7 @@ TEST_P(PcaTestDataVecSmallD, Result) data_back.data(), (params.n_col * params.n_col), raft::CompareApprox(params.tolerance), - handle.get_stream())); + resource::get_cuda_stream(handle))); } // FIXME: These tests are disabled due to driver 418+ making them fail: @@ -280,7 +299,7 @@ TEST_P(PcaTestDataVecF, Result) data2_back.data(), (params.n_col2 * params.n_col2), raft::CompareApprox(params.tolerance), - handle.get_stream())); + resource::get_cuda_stream(handle))); } typedef PcaTest PcaTestDataVecD; @@ -290,7 +309,7 @@ TEST_P(PcaTestDataVecD, Result) data2_back.data(), (params.n_col2 * params.n_col2), raft::CompareApprox(params.tolerance), - handle.get_stream())); + resource::get_cuda_stream(handle))); } INSTANTIATE_TEST_CASE_P(PcaTests, PcaTestValF, ::testing::ValuesIn(inputsf2)); diff --git a/cpp/tests/linalg/tsvd.cu b/cpp/tests/linalg/tsvd.cu index 6892fade3c..bfc211df66 100644 --- a/cpp/tests/linalg/tsvd.cu +++ b/cpp/tests/linalg/tsvd.cu @@ -5,7 +5,8 @@ #include "../test_utils.cuh" -#include +#include +#include #include #include #include @@ -40,7 +41,7 @@ class TsvdTest : public ::testing::TestWithParam> { public: TsvdTest() : params(::testing::TestWithParam>::GetParam()), - stream(handle.get_stream()), + stream(resource::get_cuda_stream(handle)), components(0, stream), components_ref(0, stream), data2(0, stream), @@ -82,7 +83,14 @@ class TsvdTest : public ::testing::TestWithParam> { else prms.algorithm = solver::COV_EIG_JACOBI; - tsvdFit(handle, data.data(), components.data(), singular_vals.data(), prms, stream); + auto input_view = raft::make_device_matrix_view( + data.data(), prms.n_rows, prms.n_cols); + auto components_view = raft::make_device_matrix_view( + components.data(), prms.n_components, prms.n_cols); + auto singular_vals_view = + raft::make_device_vector_view(singular_vals.data(), prms.n_components); + + tsvd_fit(handle, prms, input_view, components_view, singular_vals_view); } void advancedTest() @@ -122,24 +130,36 @@ class TsvdTest : public ::testing::TestWithParam> { rmm::device_uvector explained_var_ratio2(prms.n_components, stream); rmm::device_uvector singular_vals2(prms.n_components, stream); - tsvdFitTransform(handle, - data2.data(), - data2_trans.data(), - components2.data(), - explained_vars2.data(), - explained_var_ratio2.data(), - singular_vals2.data(), - prms, - stream); + auto input_view = raft::make_device_matrix_view( + data2.data(), prms.n_rows, prms.n_cols); + auto trans_view = raft::make_device_matrix_view( + data2_trans.data(), prms.n_rows, prms.n_components); + auto comp_view = raft::make_device_matrix_view( + components2.data(), prms.n_components, prms.n_cols); + auto ev_view = + raft::make_device_vector_view(explained_vars2.data(), prms.n_components); + auto evr_view = + raft::make_device_vector_view(explained_var_ratio2.data(), prms.n_components); + auto sv_view = + raft::make_device_vector_view(singular_vals2.data(), prms.n_components); + + tsvd_fit_transform(handle, prms, input_view, trans_view, comp_view, ev_view, evr_view, sv_view); data2_back.resize(len, stream); - tsvdInverseTransform( - handle, data2_trans.data(), components2.data(), data2_back.data(), prms, stream); + + auto trans_in_view = raft::make_device_matrix_view( + data2_trans.data(), prms.n_rows, prms.n_components); + auto comp_in_view = raft::make_device_matrix_view( + components2.data(), prms.n_components, prms.n_cols); + auto output_view = raft::make_device_matrix_view( + data2_back.data(), prms.n_rows, prms.n_cols); + + tsvd_inverse_transform(handle, prms, trans_in_view, comp_in_view, output_view); } protected: - raft::handle_t handle; - cudaStream_t stream = 0; + raft::device_resources handle; + cudaStream_t stream; TsvdInputs params; rmm::device_uvector components, components_ref, data2, data2_back; @@ -162,7 +182,7 @@ TEST_P(TsvdTestLeftVecF, Result) components_ref.data(), (params.n_col * params.n_col), raft::CompareApprox(params.tolerance), - handle.get_stream())); + resource::get_cuda_stream(handle))); } typedef TsvdTest TsvdTestLeftVecD; @@ -172,7 +192,7 @@ TEST_P(TsvdTestLeftVecD, Result) components_ref.data(), (params.n_col * params.n_col), raft::CompareApprox(params.tolerance), - handle.get_stream())); + resource::get_cuda_stream(handle))); } typedef TsvdTest TsvdTestDataVecF; @@ -182,7 +202,7 @@ TEST_P(TsvdTestDataVecF, Result) data2_back.data(), (params.n_col2 * params.n_col2), raft::CompareApprox(params.tolerance), - handle.get_stream())); + resource::get_cuda_stream(handle))); } typedef TsvdTest TsvdTestDataVecD; @@ -192,7 +212,7 @@ TEST_P(TsvdTestDataVecD, Result) data2_back.data(), (params.n_col2 * params.n_col2), raft::CompareApprox(params.tolerance), - handle.get_stream())); + resource::get_cuda_stream(handle))); } INSTANTIATE_TEST_CASE_P(TsvdTests, TsvdTestLeftVecF, ::testing::ValuesIn(inputsf2)); From 0c86857db3a25c56833c4281dc83a84d2b2bd597 Mon Sep 17 00:00:00 2001 From: aamijar Date: Mon, 16 Feb 2026 01:10:23 +0000 Subject: [PATCH 03/12] remove default template type --- cpp/include/raft/linalg/pca.cuh | 12 ++++++------ cpp/include/raft/linalg/tsvd.cuh | 8 ++++---- cpp/tests/linalg/pca.cu | 4 ++-- 3 files changed, 12 insertions(+), 12 deletions(-) diff --git a/cpp/include/raft/linalg/pca.cuh b/cpp/include/raft/linalg/pca.cuh index 9c305c75ff..280e833978 100644 --- a/cpp/include/raft/linalg/pca.cuh +++ b/cpp/include/raft/linalg/pca.cuh @@ -36,7 +36,7 @@ namespace raft::linalg { * @param[out] noise_vars variance of the noise. Scalar. * @param[in] flip_signs_based_on_U whether to determine signs by U (true) or V.T (false) */ -template +template void pca_fit(raft::resources const& handle, const paramsPCA& prms, raft::device_matrix_view input, @@ -45,7 +45,7 @@ void pca_fit(raft::resources const& handle, raft::device_vector_view explained_var_ratio, raft::device_vector_view singular_vals, raft::device_vector_view mu, - raft::device_scalar_view noise_vars, + raft::device_scalar_view noise_vars, bool flip_signs_based_on_U = false) { auto stream = resource::get_cuda_stream(handle); @@ -88,7 +88,7 @@ void pca_fit(raft::resources const& handle, * @param[out] noise_vars variance of the noise. Scalar. * @param[in] flip_signs_based_on_U whether to determine signs by U (true) or V.T (false) */ -template +template void pca_fit_transform(raft::resources const& handle, const paramsPCA& prms, raft::device_matrix_view input, @@ -98,7 +98,7 @@ void pca_fit_transform(raft::resources const& handle, raft::device_vector_view explained_var_ratio, raft::device_vector_view singular_vals, raft::device_vector_view mu, - raft::device_scalar_view noise_vars, + raft::device_scalar_view noise_vars, bool flip_signs_based_on_U = false) { auto stream = resource::get_cuda_stream(handle); @@ -135,7 +135,7 @@ void pca_fit_transform(raft::resources const& handle, * @param[in] mu mean of features (every column). Size n_cols. * @param[out] output the reconstructed data. Size n_rows x n_cols (col-major). */ -template +template void pca_inverse_transform(raft::resources const& handle, const paramsPCA& prms, raft::device_matrix_view trans_input, @@ -174,7 +174,7 @@ void pca_inverse_transform(raft::resources const& handle, * @param[in] mu mean value of the input data. Size n_cols. * @param[out] trans_input the transformed data. Size n_rows x n_components (col-major). */ -template +template void pca_transform(raft::resources const& handle, const paramsPCA& prms, raft::device_matrix_view input, diff --git a/cpp/include/raft/linalg/tsvd.cuh b/cpp/include/raft/linalg/tsvd.cuh index 8a54e4cf1b..9443d1c771 100644 --- a/cpp/include/raft/linalg/tsvd.cuh +++ b/cpp/include/raft/linalg/tsvd.cuh @@ -29,7 +29,7 @@ namespace raft::linalg { * @param[out] singular_vals singular values of the data. Size n_components. * @param[in] flip_signs_based_on_U whether to determine signs by U (true) or V.T (false) */ -template +template void tsvd_fit(raft::resources const& handle, const paramsTSVD& prms, raft::device_matrix_view input, @@ -70,7 +70,7 @@ void tsvd_fit(raft::resources const& handle, * @param[out] singular_vals singular values of the data. Size n_components. * @param[in] flip_signs_based_on_U whether to determine signs by U (true) or V.T (false) */ -template +template void tsvd_fit_transform(raft::resources const& handle, const paramsTSVD& prms, raft::device_matrix_view input, @@ -111,7 +111,7 @@ void tsvd_fit_transform(raft::resources const& handle, * @param[out] trans_input output that is transformed version of input. Size n_rows x n_components * (col-major). */ -template +template void tsvd_transform(raft::resources const& handle, const paramsTSVD& prms, raft::device_matrix_view input, @@ -144,7 +144,7 @@ void tsvd_transform(raft::resources const& handle, * (col-major). * @param[out] output the reconstructed data. Size n_rows x n_cols (col-major). */ -template +template void tsvd_inverse_transform(raft::resources const& handle, const paramsTSVD& prms, raft::device_matrix_view trans_input, diff --git a/cpp/tests/linalg/pca.cu b/cpp/tests/linalg/pca.cu index a0d87995fa..c0d5fe5c4f 100644 --- a/cpp/tests/linalg/pca.cu +++ b/cpp/tests/linalg/pca.cu @@ -109,7 +109,7 @@ class PcaTest : public ::testing::TestWithParam> { auto singular_vals_view = raft::make_device_vector_view(singular_vals.data(), prms.n_components); auto mu_view = raft::make_device_vector_view(mean.data(), prms.n_cols); - auto noise_vars_view = raft::make_device_scalar_view(noise_vars.data()); + auto noise_vars_view = raft::make_device_scalar_view(noise_vars.data()); pca_fit(handle, prms, @@ -173,7 +173,7 @@ class PcaTest : public ::testing::TestWithParam> { auto sv_view = raft::make_device_vector_view(singular_vals2.data(), prms.n_components); auto mu_view = raft::make_device_vector_view(mean2.data(), prms.n_cols); - auto noise_view = raft::make_device_scalar_view(noise_vars2.data()); + auto noise_view = raft::make_device_scalar_view(noise_vars2.data()); pca_fit_transform(handle, prms, From 9e285e699d687e7c5f1f60c2c122b07fae3d7ed4 Mon Sep 17 00:00:00 2001 From: aamijar Date: Mon, 16 Feb 2026 01:14:52 +0000 Subject: [PATCH 04/12] update docstring --- cpp/include/raft/linalg/detail/pca.cuh | 8 ++++---- cpp/include/raft/linalg/detail/tsvd.cuh | 10 +++++----- 2 files changed, 9 insertions(+), 9 deletions(-) diff --git a/cpp/include/raft/linalg/detail/pca.cuh b/cpp/include/raft/linalg/detail/pca.cuh index eff409c730..20ed06f59b 100644 --- a/cpp/include/raft/linalg/detail/pca.cuh +++ b/cpp/include/raft/linalg/detail/pca.cuh @@ -86,7 +86,7 @@ void truncCompExpVars(raft::resources const& handle, /** * @brief perform fit operation for the pca. Generates eigenvectors, explained vars, singular vals, * etc. - * @param[in] handle: cuml handle object + * @param[in] handle: raft::resources * @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. @@ -158,7 +158,7 @@ void pcaFit(raft::resources const& handle, /** * @brief performs transform operation for the pca. Transforms the data to eigenspace. - * @param[in] handle: the internal cuml handle object + * @param[in] handle: raft::resources * @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. @@ -208,7 +208,7 @@ void pcaTransform(raft::resources const& handle, /** * @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] handle: raft::resources * @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. @@ -259,7 +259,7 @@ void pcaInverseTransform(raft::resources const& handle, /** * @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] handle: raft::resources * @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. diff --git a/cpp/include/raft/linalg/detail/tsvd.cuh b/cpp/include/raft/linalg/detail/tsvd.cuh index 8dd95cee1c..84f1965d11 100644 --- a/cpp/include/raft/linalg/detail/tsvd.cuh +++ b/cpp/include/raft/linalg/detail/tsvd.cuh @@ -135,7 +135,7 @@ void calEig(raft::resources const& handle, /** * @defgroup sign flip for PCA and tSVD. This is used to stabilize the sign of column major eigen * vectors - * @param handle: resource handle + * @param handle: raft::resources * @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 @@ -290,7 +290,7 @@ void signFlip(math_t* input, /** * @brief perform fit operation for the tsvd. Generates eigenvectors, explained vars, singular vals, * etc. - * @param[in] handle: the internal cuml handle object + * @param[in] handle: raft::resources * @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. @@ -371,7 +371,7 @@ void tsvdFit(raft::resources const& handle, /** * @brief performs transform operation for the tsvd. Transforms the data to eigenspace. - * @param[in] handle the internal cuml handle object + * @param[in] handle raft::resources * @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 @@ -411,7 +411,7 @@ void tsvdTransform(raft::resources const& handle, /** * @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] handle raft::resources * @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. @@ -453,7 +453,7 @@ void tsvdInverseTransform(raft::resources const& handle, /** * @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] handle: raft::resources * @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. From 529514e6cfc24d646619f23b7b831f65b56316c0 Mon Sep 17 00:00:00 2001 From: aamijar Date: Mon, 16 Feb 2026 01:16:45 +0000 Subject: [PATCH 05/12] remove fixme comment --- cpp/tests/linalg/pca.cu | 2 -- 1 file changed, 2 deletions(-) diff --git a/cpp/tests/linalg/pca.cu b/cpp/tests/linalg/pca.cu index c0d5fe5c4f..beb6220432 100644 --- a/cpp/tests/linalg/pca.cu +++ b/cpp/tests/linalg/pca.cu @@ -290,8 +290,6 @@ TEST_P(PcaTestDataVecSmallD, Result) resource::get_cuda_stream(handle))); } -// FIXME: These tests are disabled due to driver 418+ making them fail: -// https://github.com/rapidsai/cuml/issues/379 typedef PcaTest PcaTestDataVecF; TEST_P(PcaTestDataVecF, Result) { From e9f6e2cd0185b9fa20b502b8fe87ff1c60e4779a Mon Sep 17 00:00:00 2001 From: aamijar Date: Mon, 16 Feb 2026 07:31:01 +0000 Subject: [PATCH 06/12] expose more tsvd functions --- cpp/include/raft/linalg/pca.cuh | 39 +++++++++++++++++++ cpp/include/raft/linalg/tsvd.cuh | 66 ++++++++++++++++++++++++++++++++ 2 files changed, 105 insertions(+) diff --git a/cpp/include/raft/linalg/pca.cuh b/cpp/include/raft/linalg/pca.cuh index 280e833978..fa721ab782 100644 --- a/cpp/include/raft/linalg/pca.cuh +++ b/cpp/include/raft/linalg/pca.cuh @@ -199,6 +199,45 @@ void pca_transform(raft::resources const& handle, stream); } +/** + * @brief Compute truncated components, explained variances, explained variance ratios, + * and noise variance from a covariance matrix. + * @tparam math_t data-type upon which the math operation will be performed + * @tparam idx_t integer type used for indexing + * @tparam enum_solver solver enum type + * @param[in] handle raft::resources + * @param[in] prms tSVD parameters (controls n_components, algorithm) + * @param[inout] in covariance matrix [n_cols x n_cols] (col-major). Overwritten. + * @param[out] components truncated eigenvectors [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] noise_vars noise variance scalar + */ +template +void trunc_comp_exp_vars(raft::resources const& handle, + const paramsTSVDTemplate& prms, + raft::device_matrix_view in, + raft::device_matrix_view components, + raft::device_vector_view explained_var, + raft::device_vector_view explained_var_ratio, + raft::device_scalar_view noise_vars) +{ + auto stream = resource::get_cuda_stream(handle); + + paramsTSVDTemplate prms_with_dims = prms; + prms_with_dims.n_rows = static_cast(in.extent(0)); + prms_with_dims.n_cols = static_cast(in.extent(1)); + + detail::truncCompExpVars(handle, + in.data_handle(), + components.data_handle(), + explained_var.data_handle(), + explained_var_ratio.data_handle(), + noise_vars.data_handle(), + prms_with_dims, + stream); +} + /** @} */ // end group pca }; // end namespace raft::linalg diff --git a/cpp/include/raft/linalg/tsvd.cuh b/cpp/include/raft/linalg/tsvd.cuh index 9443d1c771..7e6bd4731a 100644 --- a/cpp/include/raft/linalg/tsvd.cuh +++ b/cpp/include/raft/linalg/tsvd.cuh @@ -165,6 +165,72 @@ void tsvd_inverse_transform(raft::resources const& handle, stream); } +/** + * @brief Eigendecomposition helper for tSVD/PCA. Computes eigenvectors and eigenvalues + * of a symmetric matrix using either divide-and-conquer or Jacobi method. + * @tparam math_t data-type upon which the math operation will be performed + * @tparam idx_t integer type used for indexing + * @tparam enum_solver solver enum type + * @param[in] handle raft::resources + * @param[in] prms tSVD parameters (controls algorithm, tolerance, iterations) + * @param[inout] in symmetric input matrix [n_cols x n_cols] (col-major). Overwritten. + * @param[out] components eigenvectors [n_cols x n_cols] (col-major) + * @param[out] explained_var eigenvalues [n_cols] + */ +template +void cal_eig(raft::resources const& handle, + const paramsTSVDTemplate& prms, + raft::device_matrix_view in, + raft::device_matrix_view components, + raft::device_vector_view explained_var) +{ + auto stream = resource::get_cuda_stream(handle); + + paramsTSVDTemplate prms_with_dims = prms; + prms_with_dims.n_rows = static_cast(in.extent(0)); + prms_with_dims.n_cols = static_cast(in.extent(1)); + + detail::calEig(handle, + in.data_handle(), + components.data_handle(), + explained_var.data_handle(), + prms_with_dims, + stream); +} + +/** + * @brief Sign flip for PCA and tSVD. Stabilizes the sign of column-major eigenvectors. + * @tparam math_t data-type upon which the math operation will be performed + * @tparam idx_t integer type used for indexing + * @param[in] handle raft::resources + * @param[in] input input data matrix [n_samples x n_features] (col-major) + * @param[inout] components components matrix [n_components x n_features] (col-major) + * @param[in] center whether to mean-center input before computing signs + * @param[in] flip_signs_based_on_U whether to determine signs by U (true) or V.T (false) + */ +template +void sign_flip_components(raft::resources const& handle, + raft::device_matrix_view input, + raft::device_matrix_view components, + bool center, + bool flip_signs_based_on_U = false) +{ + auto stream = resource::get_cuda_stream(handle); + auto n_samples = static_cast(input.extent(0)); + auto n_features = static_cast(input.extent(1)); + auto n_components = static_cast(components.extent(0)); + + detail::signFlipComponents(handle, + input.data_handle(), + components.data_handle(), + n_samples, + n_features, + n_components, + stream, + center, + flip_signs_based_on_U); +} + /** @} */ // end group tsvd }; // end namespace raft::linalg From b85cf2533107c5b33c6b8f20da256b851cc45b28 Mon Sep 17 00:00:00 2001 From: aamijar Date: Wed, 18 Feb 2026 08:41:47 +0000 Subject: [PATCH 07/12] simplify paramsPCA --- cpp/include/raft/linalg/detail/pca.cuh | 7 ++- cpp/include/raft/linalg/detail/tsvd.cuh | 6 +-- cpp/include/raft/linalg/pca.cuh | 11 ++-- cpp/include/raft/linalg/pca_types.hpp | 70 ++++++------------------- cpp/include/raft/linalg/tsvd.cuh | 11 ++-- 5 files changed, 32 insertions(+), 73 deletions(-) diff --git a/cpp/include/raft/linalg/detail/pca.cuh b/cpp/include/raft/linalg/detail/pca.cuh index 20ed06f59b..9e7b942d39 100644 --- a/cpp/include/raft/linalg/detail/pca.cuh +++ b/cpp/include/raft/linalg/detail/pca.cuh @@ -25,14 +25,14 @@ namespace raft::linalg::detail { -template +template void truncCompExpVars(raft::resources const& handle, math_t* in, math_t* components, math_t* explained_var, math_t* explained_var_ratio, math_t* noise_vars, - const paramsTSVDTemplate& prms, + const paramsTSVD& prms, cudaStream_t stream) { auto len = prms.n_cols * prms.n_cols; @@ -40,8 +40,7 @@ void truncCompExpVars(raft::resources const& handle, rmm::device_uvector explained_var_all(prms.n_cols, stream); rmm::device_uvector explained_var_ratio_all(prms.n_cols, stream); - detail::calEig( - handle, in, components_all.data(), explained_var_all.data(), prms, stream); + detail::calEig(handle, in, components_all.data(), explained_var_all.data(), prms, stream); raft::matrix::trunc_zero_origin( handle, raft::make_device_matrix_view( diff --git a/cpp/include/raft/linalg/detail/tsvd.cuh b/cpp/include/raft/linalg/detail/tsvd.cuh index 84f1965d11..c049fc6b32 100644 --- a/cpp/include/raft/linalg/detail/tsvd.cuh +++ b/cpp/include/raft/linalg/detail/tsvd.cuh @@ -96,17 +96,17 @@ void calCompExpVarsSvd(raft::resources const& handle, raft::matrix::ratio(handle, explained_vars, explained_var_ratio, prms.n_components, stream); } -template +template void calEig(raft::resources const& handle, math_t* in, math_t* components, math_t* explained_var, - const paramsTSVDTemplate& prms, + const paramsTSVD& prms, cudaStream_t stream) { auto cusolver_handle = raft::resource::get_cusolver_dn_handle(handle); - if (prms.algorithm == enum_solver::COV_EIG_JACOBI) { + if (prms.algorithm == solver::COV_EIG_JACOBI) { raft::linalg::eigJacobi(handle, in, prms.n_cols, diff --git a/cpp/include/raft/linalg/pca.cuh b/cpp/include/raft/linalg/pca.cuh index fa721ab782..0e88837590 100644 --- a/cpp/include/raft/linalg/pca.cuh +++ b/cpp/include/raft/linalg/pca.cuh @@ -204,7 +204,6 @@ void pca_transform(raft::resources const& handle, * and noise variance from a covariance matrix. * @tparam math_t data-type upon which the math operation will be performed * @tparam idx_t integer type used for indexing - * @tparam enum_solver solver enum type * @param[in] handle raft::resources * @param[in] prms tSVD parameters (controls n_components, algorithm) * @param[inout] in covariance matrix [n_cols x n_cols] (col-major). Overwritten. @@ -213,9 +212,9 @@ void pca_transform(raft::resources const& handle, * @param[out] explained_var_ratio explained variance ratios [n_components] * @param[out] noise_vars noise variance scalar */ -template +template void trunc_comp_exp_vars(raft::resources const& handle, - const paramsTSVDTemplate& prms, + const paramsTSVD& prms, raft::device_matrix_view in, raft::device_matrix_view components, raft::device_vector_view explained_var, @@ -224,9 +223,9 @@ void trunc_comp_exp_vars(raft::resources const& handle, { auto stream = resource::get_cuda_stream(handle); - paramsTSVDTemplate prms_with_dims = prms; - prms_with_dims.n_rows = static_cast(in.extent(0)); - prms_with_dims.n_cols = static_cast(in.extent(1)); + paramsTSVD prms_with_dims = prms; + prms_with_dims.n_rows = static_cast(in.extent(0)); + prms_with_dims.n_cols = static_cast(in.extent(1)); detail::truncCompExpVars(handle, in.data_handle(), diff --git a/cpp/include/raft/linalg/pca_types.hpp b/cpp/include/raft/linalg/pca_types.hpp index ec70b1a267..1eba5e7fdb 100644 --- a/cpp/include/raft/linalg/pca_types.hpp +++ b/cpp/include/raft/linalg/pca_types.hpp @@ -10,69 +10,31 @@ namespace raft::linalg { /** - * @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 + * @brief Solver algorithm for PCA/TSVD eigen decomposition. + * + * @param COV_EIG_DQ covariance + divide-and-conquer eigen decomposition for symmetric matrices + * @param COV_EIG_JACOBI covariance + Jacobi eigen decomposition for symmetric matrices */ enum class solver : int { COV_EIG_DQ, COV_EIG_JACOBI, }; -class params { - public: - std::size_t n_rows; - std::size_t n_cols; - int gpu_id = 0; +/** @brief Parameters for TSVD (and base for PCA). */ +struct paramsTSVD { + std::size_t n_rows = 0; + std::size_t n_cols = 0; + int gpu_id = 0; + float tol = 0.0; + uint64_t n_iterations = 15; + uint64_t n_components = 1; + solver algorithm = solver::COV_EIG_DQ; }; -class paramsSolver : public params { - public: - // math_t tol = 0.0; - float tol = 0.0; - std::uint32_t n_iterations = 15; - int verbose = 0; -}; - -template -class paramsTSVDTemplate : public paramsSolver { - public: - std::size_t n_components = 1; - enum_solver algorithm = enum_solver::COV_EIG_DQ; -}; - -/** - * @brief structure for pca parameters. Ref: - * http://scikit-learn.org/stable/modules/generated/sklearn.decomposition.PCA.html - * @param n_components: Number of components to keep. if n_components is not set all components are - * kept: - * @param copy: 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. - * @param whiten: 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. - * @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 verbose: 0: no error message printing, 1: print error messages - */ - -template -class paramsPCATemplate : public paramsTSVDTemplate { - public: - bool copy = true; // TODO unused, see #2830 and #2833 +/** @brief Parameters for PCA (extends TSVD with whitening / copy controls). */ +struct paramsPCA : paramsTSVD { + bool copy = true; bool whiten = false; }; -typedef paramsTSVDTemplate<> paramsTSVD; -typedef paramsPCATemplate<> paramsPCA; - -enum class mg_solver { COV_EIG_DQ, COV_EIG_JACOBI }; - -typedef paramsPCATemplate paramsPCAMG; -typedef paramsTSVDTemplate paramsTSVDMG; - }; // end namespace raft::linalg diff --git a/cpp/include/raft/linalg/tsvd.cuh b/cpp/include/raft/linalg/tsvd.cuh index 7e6bd4731a..fe4eb9437b 100644 --- a/cpp/include/raft/linalg/tsvd.cuh +++ b/cpp/include/raft/linalg/tsvd.cuh @@ -170,25 +170,24 @@ void tsvd_inverse_transform(raft::resources const& handle, * of a symmetric matrix using either divide-and-conquer or Jacobi method. * @tparam math_t data-type upon which the math operation will be performed * @tparam idx_t integer type used for indexing - * @tparam enum_solver solver enum type * @param[in] handle raft::resources * @param[in] prms tSVD parameters (controls algorithm, tolerance, iterations) * @param[inout] in symmetric input matrix [n_cols x n_cols] (col-major). Overwritten. * @param[out] components eigenvectors [n_cols x n_cols] (col-major) * @param[out] explained_var eigenvalues [n_cols] */ -template +template void cal_eig(raft::resources const& handle, - const paramsTSVDTemplate& prms, + const paramsTSVD& prms, raft::device_matrix_view in, raft::device_matrix_view components, raft::device_vector_view explained_var) { auto stream = resource::get_cuda_stream(handle); - paramsTSVDTemplate prms_with_dims = prms; - prms_with_dims.n_rows = static_cast(in.extent(0)); - prms_with_dims.n_cols = static_cast(in.extent(1)); + paramsTSVD prms_with_dims = prms; + prms_with_dims.n_rows = static_cast(in.extent(0)); + prms_with_dims.n_cols = static_cast(in.extent(1)); detail::calEig(handle, in.data_handle(), From 5805ead5d034ce4294ecfc230cce7009bc04e9f9 Mon Sep 17 00:00:00 2001 From: Anupam <54245698+aamijar@users.noreply.github.com> Date: Tue, 10 Mar 2026 17:44:51 -0700 Subject: [PATCH 08/12] Update cpp/include/raft/linalg/pca.cuh Co-authored-by: Jinsol Park --- cpp/include/raft/linalg/pca.cuh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cpp/include/raft/linalg/pca.cuh b/cpp/include/raft/linalg/pca.cuh index 0e88837590..8766aa8528 100644 --- a/cpp/include/raft/linalg/pca.cuh +++ b/cpp/include/raft/linalg/pca.cuh @@ -21,7 +21,7 @@ namespace raft::linalg { * @brief perform fit operation for PCA. Generates eigenvectors, explained vars, singular vals, etc. * @tparam math_t data-type upon which the math operation will be performed * @tparam idx_t integer type used for indexing - * @param[in] handle raft::resources + * @param[in] handle: raft::resources * @param[in] prms PCA parameters (n_components, algorithm, whiten, etc.) * @param[inout] input the data is fitted to PCA. Size n_rows x n_cols (col-major). Modified * temporarily during computation. From 5945c1f3cac51d8e72f123f443a7f2ccb80297b6 Mon Sep 17 00:00:00 2001 From: aamijar Date: Wed, 18 Mar 2026 02:12:40 +0000 Subject: [PATCH 09/12] clean up detail apis --- cpp/include/raft/linalg/detail/pca.cuh | 372 ++++++++++--------- cpp/include/raft/linalg/detail/tsvd.cuh | 471 +++++++++++++----------- cpp/include/raft/linalg/pca.cuh | 97 ++--- cpp/include/raft/linalg/tsvd.cuh | 90 +---- 4 files changed, 478 insertions(+), 552 deletions(-) diff --git a/cpp/include/raft/linalg/detail/pca.cuh b/cpp/include/raft/linalg/detail/pca.cuh index 9e7b942d39..6a8ed8c17d 100644 --- a/cpp/include/raft/linalg/detail/pca.cuh +++ b/cpp/include/raft/linalg/detail/pca.cuh @@ -5,7 +5,9 @@ #pragma once +#include #include +#include #include #include #include @@ -25,59 +27,68 @@ namespace raft::linalg::detail { -template +template void truncCompExpVars(raft::resources const& handle, - math_t* in, - math_t* components, - math_t* explained_var, - math_t* explained_var_ratio, - math_t* noise_vars, - const paramsTSVD& prms, - cudaStream_t stream) + raft::device_matrix_view in, + raft::device_matrix_view components, + raft::device_vector_view explained_var, + raft::device_vector_view explained_var_ratio, + raft::device_scalar_view noise_vars, + const paramsTSVD& prms) { - auto len = prms.n_cols * prms.n_cols; + auto stream = resource::get_cuda_stream(handle); + + auto n_cols = in.extent(0); + auto n_components = components.extent(0); + + auto len = static_cast(n_cols * 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); + rmm::device_uvector explained_var_all(static_cast(n_cols), stream); + rmm::device_uvector explained_var_ratio_all(static_cast(n_cols), stream); - detail::calEig(handle, in, components_all.data(), explained_var_all.data(), prms, stream); + detail::calEig( + handle, + in, + raft::make_device_matrix_view( + components_all.data(), n_cols, n_cols), + raft::make_device_vector_view(explained_var_all.data(), n_cols), + prms); 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::make_device_matrix_view( + components_all.data(), n_cols, n_cols), + raft::make_device_matrix_view( + components.data_handle(), n_components, 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::make_device_matrix_view( + explained_var_all.data(), n_cols, idx_t(1)), + raft::make_device_matrix_view( + explained_var_ratio_all.data(), n_cols, idx_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::make_device_matrix_view( + explained_var_all.data(), n_cols, idx_t(1)), + raft::make_device_matrix_view( + explained_var.data_handle(), n_components, idx_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, + raft::make_device_matrix_view( + explained_var_ratio_all.data(), n_cols, idx_t(1)), + raft::make_device_matrix_view( + explained_var_ratio.data_handle(), n_components, idx_t(1))); + + if (static_cast(n_components) < static_cast(n_cols) && + static_cast(n_components) < prms.n_rows) { + raft::stats::mean(noise_vars.data_handle(), + explained_var_all.data() + static_cast(n_components), std::size_t{1}, - prms.n_cols - prms.n_components, + static_cast(n_cols - n_components), false, stream); } else { raft::matrix::fill( handle, - raft::make_device_vector_view(noise_vars, std::size_t(1)), + raft::make_device_vector_view(noise_vars.data_handle(), idx_t(1)), math_t{0}); } } @@ -86,205 +97,219 @@ void truncCompExpVars(raft::resources const& handle, * @brief perform fit operation for the pca. Generates eigenvectors, explained vars, singular vals, * etc. * @param[in] handle: raft::resources - * @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[inout] input: the data is fitted to PCA. Size n_rows x n_cols (col-major). + * @param[out] components: the principal components. Size n_components x n_cols (col-major). + * @param[out] explained_var: explained variances. Size n_components. + * @param[out] explained_var_ratio: ratio of explained to total variance. Size n_components. + * @param[out] singular_vals: singular values. Size n_components. + * @param[out] mu: mean of all features. Size n_cols. + * @param[out] noise_vars: noise variance scalar. * @param[in] prms: data structure that includes all the parameters from input size to algorithm. - * @param[in] stream cuda stream + * @param[in] flip_signs_based_on_U whether to determine signs by U (true) or V.T (false) */ -template +template void pcaFit(raft::resources const& 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, + 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, const paramsPCA& prms, - cudaStream_t stream, bool flip_signs_based_on_U = false) { + auto stream = resource::get_cuda_stream(handle); auto cublas_handle = raft::resource::get_cublas_handle(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"); + auto n_rows = input.extent(0); + auto n_cols = input.extent(1); + + ASSERT(n_cols > 1, "Parameter n_cols: number of columns cannot be less than two"); + ASSERT(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; + auto n_components = static_cast(prms.n_components); + if (n_components > n_cols) n_components = n_cols; - raft::stats::mean(mu, input, prms.n_cols, prms.n_rows, false, stream); + raft::stats::mean(mu.data_handle(), input.data_handle(), n_cols, n_rows, false, stream); - auto len = prms.n_cols * prms.n_cols; + auto len = static_cast(n_cols * 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); - detail::truncCompExpVars( - handle, cov.data(), components, explained_var, explained_var_ratio, noise_vars, prms, stream); + handle, cov.data(), input.data_handle(), mu.data_handle(), n_cols, n_rows, true, true, stream); + + paramsPCA prms_with_rows = prms; + prms_with_rows.n_rows = static_cast(n_rows); + prms_with_rows.n_cols = static_cast(n_cols); - math_t scalar = (prms.n_rows - 1); - raft::matrix::weighted_sqrt( + detail::truncCompExpVars( 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); - - detail::signFlipComponents(handle, - input, - components, - prms.n_rows, - prms.n_cols, - prms.n_components, - stream, - true, - flip_signs_based_on_U); + raft::make_device_matrix_view(cov.data(), n_cols, n_cols), + components, + explained_var, + explained_var_ratio, + noise_vars, + prms_with_rows); + + math_t scalar = (n_rows - 1); + raft::matrix::weighted_sqrt(handle, + raft::make_device_matrix_view( + explained_var.data_handle(), idx_t(1), n_components), + raft::make_device_matrix_view( + singular_vals.data_handle(), idx_t(1), n_components), + raft::make_host_scalar_view(&scalar), + true); + + raft::stats::meanAdd( + input.data_handle(), input.data_handle(), mu.data_handle(), n_cols, n_rows, stream); + + detail::signFlipComponents(handle, input, components, true, flip_signs_based_on_U); } /** * @brief performs transform operation for the pca. Transforms the data to eigenspace. * @param[in] handle: raft::resources - * @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[inout] input: the data to transform. Size n_rows x n_cols (col-major). Modified + * temporarily (mean-centered then restored). + * @param[in] components: principal components. Size n_components x n_cols (col-major). + * @param[out] trans_input: the transformed data. Size n_rows x n_components (col-major). + * @param[in] singular_vals: singular values. Size n_components. + * @param[in] mu: mean of features. Size n_cols. * @param[in] prms: data structure that includes all the parameters from input size to algorithm. - * @param[in] stream cuda stream */ -template +template void pcaTransform(raft::resources const& handle, - math_t* input, - math_t* components, - math_t* trans_input, - math_t* singular_vals, - math_t* mu, - const paramsPCA& prms, - cudaStream_t stream) + raft::device_matrix_view input, + raft::device_matrix_view components, + raft::device_matrix_view trans_input, + raft::device_vector_view singular_vals, + raft::device_vector_view mu, + const paramsPCA& prms) { - 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 stream = resource::get_cuda_stream(handle); + + auto n_rows = input.extent(0); + auto n_cols = input.extent(1); + auto n_components = components.extent(0); + + ASSERT(n_cols > 1, "Parameter n_cols: number of columns cannot be less than two"); + ASSERT(n_rows > 0, "Parameter n_rows: number of rows cannot be less than one"); + ASSERT(n_components > 0, "Parameter n_components: number of components cannot be less than one"); - auto components_len = prms.n_cols * prms.n_components; + auto components_len = static_cast(n_cols * n_components); rmm::device_uvector components_copy{components_len, stream}; - raft::copy(components_copy.data(), components, prms.n_cols * prms.n_components, stream); + raft::copy(components_copy.data(), components.data_handle(), components_len, 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); + math_t scalar = math_t(sqrt(n_rows - 1)); + raft::linalg::scalarMultiply( + components_copy.data(), components_copy.data(), scalar, components_len, 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::make_device_matrix_view( + components_copy.data(), n_cols, n_components), + raft::make_device_vector_view(singular_vals.data_handle(), + n_components)); } - raft::stats::meanCenter(input, input, mu, prms.n_cols, prms.n_rows, stream); - detail::tsvdTransform(handle, input, components_copy.data(), trans_input, prms, stream); - raft::stats::meanAdd(input, input, mu, prms.n_cols, prms.n_rows, stream); + raft::stats::meanCenter( + input.data_handle(), input.data_handle(), mu.data_handle(), n_cols, n_rows, stream); + detail::tsvdTransform(handle, + input, + raft::make_device_matrix_view( + components_copy.data(), n_components, n_cols), + trans_input, + prms); + raft::stats::meanAdd( + input.data_handle(), input.data_handle(), mu.data_handle(), n_cols, n_rows, stream); } /** * @brief performs inverse transform operation for the pca. Transforms the transformed data back to * original data. * @param[in] handle: raft::resources - * @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] trans_input: the transformed data. Size n_rows x n_components (col-major). + * @param[in] components: principal components. Size n_components x n_cols (col-major). + * @param[in] singular_vals: singular values. Size n_components. + * @param[in] mu: mean of features. Size n_cols. + * @param[out] input: the reconstructed data. Size n_rows x n_cols (col-major). * @param[in] prms: data structure that includes all the parameters from input size to algorithm. - * @param[in] stream cuda stream */ -template +template void pcaInverseTransform(raft::resources const& handle, - math_t* trans_input, - math_t* components, - math_t* singular_vals, - math_t* mu, - math_t* input, - const paramsPCA& prms, - cudaStream_t stream) + 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 input, + const paramsPCA& prms) { - 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 stream = resource::get_cuda_stream(handle); + + auto n_rows = input.extent(0); + auto n_cols = input.extent(1); + auto n_components = components.extent(0); + + ASSERT(n_cols > 1, "Parameter n_cols: number of columns cannot be less than two"); + ASSERT(n_rows > 0, "Parameter n_rows: number of rows cannot be less than one"); + ASSERT(n_components > 0, "Parameter n_components: number of components cannot be less than one"); - auto components_len = prms.n_cols * prms.n_components; + auto components_len = static_cast(n_cols * n_components); rmm::device_uvector components_copy{components_len, stream}; - raft::copy(components_copy.data(), components, prms.n_cols * prms.n_components, stream); + raft::copy(components_copy.data(), components.data_handle(), components_len, 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); + math_t sqrt_n_samples = sqrt(n_rows - 1); + math_t scalar = n_rows - 1 > 0 ? math_t(1 / sqrt_n_samples) : 0; + raft::linalg::scalarMultiply( + components_copy.data(), components_copy.data(), scalar, components_len, 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)); + raft::make_device_matrix_view( + components_copy.data(), n_cols, n_components), + raft::make_device_vector_view(singular_vals.data_handle(), + n_components)); } - detail::tsvdInverseTransform(handle, trans_input, components_copy.data(), input, prms, stream); - raft::stats::meanAdd(input, input, mu, prms.n_cols, prms.n_rows, stream); + detail::tsvdInverseTransform(handle, + trans_input, + raft::make_device_matrix_view( + components_copy.data(), n_components, n_cols), + input, + prms); + raft::stats::meanAdd( + input.data_handle(), input.data_handle(), mu.data_handle(), n_cols, n_rows, stream); } /** * @brief perform fit and transform operations for the pca. Generates transformed data, * eigenvectors, explained vars, singular vals, etc. * @param[in] handle: raft::resources - * @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[inout] input: the data is fitted to PCA. Size n_rows x n_cols (col-major). + * @param[out] trans_input: the transformed data. Size n_rows x n_components (col-major). + * @param[out] components: the principal components. Size n_components x n_cols (col-major). + * @param[out] explained_var: explained variances. Size n_components. + * @param[out] explained_var_ratio: ratio of explained to total variance. Size n_components. + * @param[out] singular_vals: singular values. Size n_components. + * @param[out] mu: mean of all features. Size n_cols. + * @param[out] noise_vars: noise variance scalar. * @param[in] prms: data structure that includes all the parameters from input size to algorithm. - * @param[in] stream cuda stream + * @param[in] flip_signs_based_on_U whether to determine signs by U (true) or V.T (false) */ -template +template void pcaFitTransform(raft::resources const& 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, + 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, const paramsPCA& prms, - cudaStream_t stream, bool flip_signs_based_on_U = false) { detail::pcaFit(handle, @@ -296,9 +321,8 @@ void pcaFitTransform(raft::resources const& handle, mu, noise_vars, prms, - stream, flip_signs_based_on_U); - detail::pcaTransform(handle, input, components, trans_input, singular_vals, mu, prms, stream); + detail::pcaTransform(handle, input, components, trans_input, singular_vals, mu, prms); } }; // end namespace raft::linalg::detail diff --git a/cpp/include/raft/linalg/detail/tsvd.cuh b/cpp/include/raft/linalg/detail/tsvd.cuh index c049fc6b32..5494190e49 100644 --- a/cpp/include/raft/linalg/detail/tsvd.cuh +++ b/cpp/include/raft/linalg/detail/tsvd.cuh @@ -5,8 +5,10 @@ #pragma once +#include #include #include +#include #include #include #include @@ -37,43 +39,47 @@ namespace raft::linalg::detail { -template +template void calCompExpVarsSvd(raft::resources const& 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) + raft::device_matrix_view in, + raft::device_matrix_view components, + raft::device_vector_view singular_vals, + raft::device_vector_view explained_vars, + raft::device_vector_view explained_var_ratio, + const paramsTSVD& prms) { + auto stream = resource::get_cuda_stream(handle); auto cusolver_handle = raft::resource::get_cusolver_dn_handle(handle); auto cublas_handle = raft::resource::get_cublas_handle(handle); - auto diff = prms.n_cols - prms.n_components; - math_t ratio = math_t(diff) / math_t(prms.n_cols); + auto n_rows = in.extent(0); + auto n_cols = in.extent(1); + auto n_components = components.extent(0); + + auto diff = n_cols - n_components; + math_t ratio = math_t(diff) / math_t(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)); + std::size_t p = static_cast(math_t(0.1) * math_t(n_cols)); 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, + auto total_random_vecs = static_cast(n_components) + p; + ASSERT(total_random_vecs < static_cast(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); + rmm::device_uvector components_temp(static_cast(n_cols * n_components), + stream); math_t* left_eigvec = nullptr; raft::linalg::rsvdFixedRank(handle, - in, - prms.n_rows, - prms.n_cols, - singular_vals, + in.data_handle(), + n_rows, + n_cols, + singular_vals.data_handle(), left_eigvec, components_temp.data(), - prms.n_components, + n_components, p, true, false, @@ -84,105 +90,103 @@ void calCompExpVarsSvd(raft::resources const& handle, 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); + handle, components_temp.data(), components.data_handle(), n_cols, n_components, stream); + + raft::matrix::weighted_power(handle, + raft::make_device_matrix_view( + singular_vals.data_handle(), idx_t(1), n_components), + raft::make_device_matrix_view( + explained_vars.data_handle(), idx_t(1), n_components), + math_t(1)); + raft::matrix::ratio( + handle, explained_vars.data_handle(), explained_var_ratio.data_handle(), n_components, stream); } -template +template void calEig(raft::resources const& handle, - math_t* in, - math_t* components, - math_t* explained_var, - const paramsTSVD& prms, - cudaStream_t stream) + raft::device_matrix_view in, + raft::device_matrix_view components, + raft::device_vector_view explained_var, + const paramsTSVD& prms) { + auto stream = resource::get_cuda_stream(handle); auto cusolver_handle = raft::resource::get_cusolver_dn_handle(handle); + auto n_cols = in.extent(0); + if (prms.algorithm == solver::COV_EIG_JACOBI) { raft::linalg::eigJacobi(handle, - in, - prms.n_cols, - prms.n_cols, - components, - explained_var, + in.data_handle(), + n_cols, + n_cols, + components.data_handle(), + explained_var.data_handle(), 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::linalg::eigDC(handle, + in.data_handle(), + n_cols, + n_cols, + components.data_handle(), + explained_var.data_handle(), + 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::make_device_matrix_view( + components.data_handle(), n_cols, n_cols)); + raft::linalg::transpose(components.data_handle(), n_cols, stream); raft::matrix::row_reverse(handle_stream_zero, - raft::make_device_matrix_view( - explained_var, prms.n_cols, std::size_t(1))); + raft::make_device_matrix_view( + explained_var.data_handle(), n_cols, idx_t(1))); } /** * @defgroup sign flip for PCA and tSVD. This is used to stabilize the sign of column major eigen * vectors * @param handle: raft::resources - * @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 input: input data [n_samples x n_features] (col-major) + * @param components: components matrix [n_components x n_features] (col-major) + * @param center whether to mean-center input before computing signs * @param flip_signs_based_on_U whether to determine signs by U (true) or V.T (false) * @{ */ -template +template void signFlipComponents(raft::resources const& 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, + raft::device_matrix_view input, + raft::device_matrix_view components, 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 + auto stream = resource::get_cuda_stream(handle); + auto n_samples = input.extent(0); + auto n_features = input.extent(1); + auto n_components = components.extent(0); + + rmm::device_uvector max_vals(static_cast(n_components), stream); + auto components_view = raft::make_device_matrix_view( + components.data_handle(), n_components, n_features); + auto max_vals_view = raft::make_device_vector_view(max_vals.data(), n_components); + 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); + rmm::device_uvector col_means(static_cast(n_features), stream); + raft::stats::mean( + col_means.data(), input.data_handle(), n_features, n_samples, stream); raft::stats::meanCenter( - input, input, col_means.data(), n_features, n_samples, stream); + input.data_handle(), input.data_handle(), col_means.data(), n_features, n_samples, stream); } - rmm::device_uvector US(n_samples * n_components, stream); + rmm::device_uvector US(static_cast(n_samples * n_components), stream); raft::linalg::gemm(handle, - input, + input.data_handle(), n_samples, n_features, - components, + components.data_handle(), US.data(), n_samples, n_components, @@ -209,7 +213,7 @@ void signFlipComponents(raft::resources const& handle, } else { raft::linalg::reduce( max_vals.data(), - components, + components.data_handle(), n_features, n_components, math_t(0), @@ -224,13 +228,12 @@ void signFlipComponents(raft::resources const& handle, 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; + auto row = idx % n_components; + auto column = idx / n_components; return (max_vals_view(row) < math_t(0)) ? (-components_view(row, column)) : components_view(row, column); }); @@ -239,34 +242,35 @@ void signFlipComponents(raft::resources const& handle, /** * @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 + * @param handle: raft::resources + * @param input: input matrix [n_rows x n_cols] (col-major). Modified in place. + * @param components: components matrix [n_rows x n_cols_comp] (col-major). Modified in place. * @{ */ -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) +template +void signFlip(raft::resources const& handle, + raft::device_matrix_view input, + raft::device_matrix_view components) { - auto counting = thrust::make_counting_iterator(0); - auto m = n_rows; + auto stream = resource::get_cuda_stream(handle); + auto n_rows = input.extent(0); + auto n_cols = input.extent(1); + auto n_cols_comp = components.extent(1); + + auto* input_ptr = input.data_handle(); + auto* components_ptr = components.data_handle(); + 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) { + rmm::exec_policy(stream), counting, counting + n_cols, [=] __device__(idx_t idx) { auto d_i = idx * m; auto end = d_i + m; - math_t max = 0.0; - std::size_t max_index = 0; + math_t max = 0.0; + idx_t max_index = 0; for (auto i = d_i; i < end; i++) { - math_t val = input[i]; + math_t val = input_ptr[i]; if (val < 0.0) { val = -val; } if (val > max) { max = val; @@ -274,14 +278,14 @@ void signFlip(math_t* input, } } - if (input[max_index] < 0.0) { + if (input_ptr[max_index] < 0.0) { for (auto i = d_i; i < end; i++) { - input[i] = -input[i]; + input_ptr[i] = -input_ptr[i]; } auto len = n_cols * n_cols_comp; for (auto i = idx; i < len; i = i + n_cols) { - components[i] = -components[i]; + components_ptr[i] = -components_ptr[i]; } } }); @@ -291,45 +295,47 @@ void signFlip(math_t* input, * @brief perform fit operation for the tsvd. Generates eigenvectors, explained vars, singular vals, * etc. * @param[in] handle: raft::resources - * @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] input: the data is fitted to tSVD. Size n_rows x n_cols (col-major). + * @param[out] components: the principal components. Size n_components x n_cols (col-major). + * @param[out] singular_vals: singular values of the data. Size n_components. * @param[in] prms: data structure that includes all the parameters from input size to algorithm. - * @param[in] stream cuda stream + * @param[in] flip_signs_based_on_U whether to determine signs by U (true) or V.T (false) */ -template +template void tsvdFit(raft::resources const& handle, - math_t* input, - math_t* components, - math_t* singular_vals, + raft::device_matrix_view input, + raft::device_matrix_view components, + raft::device_vector_view singular_vals, const paramsTSVD& prms, - cudaStream_t stream, bool flip_signs_based_on_U = false) { + auto stream = resource::get_cuda_stream(handle); auto cublas_handle = raft::resource::get_cublas_handle(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"); + auto n_rows = input.extent(0); + auto n_cols = input.extent(1); + + ASSERT(n_cols > 1, "Parameter n_cols: number of columns cannot be less than two"); + ASSERT(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; + auto n_components = static_cast(prms.n_components); + if (n_components > n_cols) n_components = n_cols; - size_t len = prms.n_cols * prms.n_cols; + auto len = static_cast(n_cols * 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.data_handle(), + n_rows, + n_cols, + input.data_handle(), input_cross_mult.data(), - prms.n_cols, - prms.n_cols, + n_cols, + n_cols, CUBLAS_OP_T, CUBLAS_OP_N, alpha, @@ -337,70 +343,74 @@ void tsvdFit(raft::resources const& handle, stream); rmm::device_uvector components_all(len, stream); - rmm::device_uvector explained_var_all(prms.n_cols, stream); + rmm::device_uvector explained_var_all(static_cast(n_cols), stream); - detail::calEig( - handle, input_cross_mult.data(), components_all.data(), explained_var_all.data(), prms, stream); + detail::calEig(handle, + raft::make_device_matrix_view( + input_cross_mult.data(), n_cols, n_cols), + raft::make_device_matrix_view( + components_all.data(), n_cols, n_cols), + raft::make_device_vector_view(explained_var_all.data(), n_cols), + prms); 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)); + raft::make_device_matrix_view( + components_all.data(), n_cols, n_cols), + raft::make_device_matrix_view( + components.data_handle(), n_components, 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); + raft::matrix::weighted_sqrt(handle, + raft::make_device_matrix_view( + explained_var_all.data(), idx_t(1), n_components), + raft::make_device_matrix_view( + singular_vals.data_handle(), idx_t(1), n_components), + raft::make_host_scalar_view(&scalar)); + + detail::signFlipComponents(handle, + input, + raft::make_device_matrix_view( + components.data_handle(), n_components, n_cols), + false, + flip_signs_based_on_U); } /** * @brief performs transform operation for the tsvd. Transforms the data to eigenspace. * @param[in] handle raft::resources - * @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] input: the data to transform. Size n_rows x n_cols (col-major). + * @param[in] components: principal components. Size n_components x n_cols (col-major). + * @param[out] trans_input: transformed output. Size n_rows x n_components (col-major). * @param[in] prms: data structure that includes all the parameters from input size to algorithm. - * @param[in] stream cuda stream */ -template +template void tsvdTransform(raft::resources const& handle, - math_t* input, - math_t* components, - math_t* trans_input, - const paramsTSVD& prms, - cudaStream_t stream) + raft::device_matrix_view input, + raft::device_matrix_view components, + raft::device_matrix_view trans_input, + const paramsTSVD& prms) { - 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 stream = resource::get_cuda_stream(handle); + + auto n_rows = input.extent(0); + auto n_cols = input.extent(1); + auto n_components = components.extent(0); + + ASSERT(n_cols > 1, "Parameter n_cols: number of columns cannot be less than two"); + ASSERT(n_rows > 0, "Parameter n_rows: number of rows cannot be less than one"); + ASSERT(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, + input.data_handle(), + n_rows, + n_cols, + components.data_handle(), + trans_input.data_handle(), + n_rows, + n_components, CUBLAS_OP_N, CUBLAS_OP_T, alpha, @@ -412,37 +422,39 @@ void tsvdTransform(raft::resources const& handle, * @brief performs inverse transform operation for the tsvd. Transforms the transformed data back to * original data. * @param[in] handle raft::resources - * @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] trans_input: the transformed data. Size n_rows x n_components (col-major). + * @param[in] components: principal components. Size n_components x n_cols (col-major). + * @param[out] input: reconstructed output. Size n_rows x n_cols (col-major). * @param[in] prms: data structure that includes all the parameters from input size to algorithm. - * @param[in] stream cuda stream */ -template +template void tsvdInverseTransform(raft::resources const& handle, - math_t* trans_input, - math_t* components, - math_t* input, - const paramsTSVD& prms, - cudaStream_t stream) + raft::device_matrix_view trans_input, + raft::device_matrix_view components, + raft::device_matrix_view input, + const paramsTSVD& prms) { - 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"); + auto stream = resource::get_cuda_stream(handle); + + auto n_rows = input.extent(0); + auto n_cols = input.extent(1); + auto n_components = components.extent(0); + + ASSERT(n_cols > 1, "Parameter n_cols: number of columns cannot be less than one"); + ASSERT(n_rows > 0, "Parameter n_rows: number of rows cannot be less than one"); + ASSERT(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, + trans_input.data_handle(), + n_rows, + n_components, + components.data_handle(), + input.data_handle(), + n_rows, + n_cols, CUBLAS_OP_N, CUBLAS_OP_N, alpha, @@ -454,47 +466,56 @@ void tsvdInverseTransform(raft::resources const& handle, * @brief performs fit and transform operations for the tsvd. Generates transformed data, * eigenvectors, explained vars, singular vals, etc. * @param[in] handle: raft::resources - * @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] input: the data is fitted to tSVD. Size n_rows x n_cols (col-major). + * @param[out] trans_input: the transformed data. Size n_rows x n_components (col-major). + * @param[out] components: the principal components. Size n_components x n_cols (col-major). + * @param[out] explained_var: explained variances. Size n_components. + * @param[out] explained_var_ratio: ratio of explained variance to total. Size n_components. + * @param[out] singular_vals: singular values of the data. Size n_components. * @param[in] prms: data structure that includes all the parameters from input size to algorithm. - * @param[in] stream cuda stream + * @param[in] flip_signs_based_on_U whether to determine signs by U (true) or V.T (false) */ -template +template void tsvdFitTransform(raft::resources const& handle, - math_t* input, - math_t* trans_input, - math_t* components, - math_t* explained_var, - math_t* explained_var_ratio, - math_t* singular_vals, + 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, const paramsTSVD& prms, - cudaStream_t stream, bool flip_signs_based_on_U = false) { - detail::tsvdFit(handle, input, components, singular_vals, prms, stream, flip_signs_based_on_U); - detail::tsvdTransform(handle, input, components, trans_input, prms, stream); + auto stream = resource::get_cuda_stream(handle); - 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); + auto n_rows = input.extent(0); + auto n_cols = input.extent(1); + auto n_components = components.extent(0); - rmm::device_uvector mu(prms.n_cols, stream); - rmm::device_uvector vars(prms.n_cols, stream); + detail::tsvdFit(handle, input, components, singular_vals, prms, flip_signs_based_on_U); + detail::tsvdTransform(handle, input, components, trans_input, prms); - 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_uvector mu_trans(static_cast(n_components), stream); + raft::stats::mean( + mu_trans.data(), trans_input.data_handle(), n_components, n_rows, false, stream); + raft::stats::vars(explained_var.data_handle(), + trans_input.data_handle(), + mu_trans.data(), + n_components, + n_rows, + false, + stream); + + rmm::device_uvector mu(static_cast(n_cols), stream); + rmm::device_uvector vars(static_cast(n_cols), stream); + + raft::stats::mean(mu.data(), input.data_handle(), n_cols, n_rows, false, stream); + raft::stats::vars( + vars.data(), input.data_handle(), mu.data(), n_cols, 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); + raft::stats::sum( + total_vars.data(), vars.data(), std::size_t(1), static_cast(n_cols), stream); math_t total_vars_h; raft::update_host(&total_vars_h, total_vars.data(), 1, stream); @@ -502,7 +523,7 @@ void tsvdFitTransform(raft::resources const& handle, math_t scalar = math_t(1) / total_vars_h; raft::linalg::scalarMultiply( - explained_var_ratio, explained_var, scalar, prms.n_components, stream); + explained_var_ratio.data_handle(), explained_var.data_handle(), scalar, n_components, stream); } }; // end namespace raft::linalg::detail diff --git a/cpp/include/raft/linalg/pca.cuh b/cpp/include/raft/linalg/pca.cuh index 8766aa8528..5136476884 100644 --- a/cpp/include/raft/linalg/pca.cuh +++ b/cpp/include/raft/linalg/pca.cuh @@ -8,7 +8,6 @@ #include "detail/pca.cuh" #include -#include namespace raft::linalg { @@ -48,22 +47,15 @@ void pca_fit(raft::resources const& handle, raft::device_scalar_view noise_vars, bool flip_signs_based_on_U = false) { - auto stream = resource::get_cuda_stream(handle); - - paramsPCA prms_with_dims = prms; - prms_with_dims.n_rows = static_cast(input.extent(0)); - prms_with_dims.n_cols = static_cast(input.extent(1)); - detail::pcaFit(handle, - input.data_handle(), - components.data_handle(), - explained_var.data_handle(), - explained_var_ratio.data_handle(), - singular_vals.data_handle(), - mu.data_handle(), - noise_vars.data_handle(), - prms_with_dims, - stream, + input, + components, + explained_var, + explained_var_ratio, + singular_vals, + mu, + noise_vars, + prms, flip_signs_based_on_U); } @@ -101,23 +93,16 @@ void pca_fit_transform(raft::resources const& handle, raft::device_scalar_view noise_vars, bool flip_signs_based_on_U = false) { - auto stream = resource::get_cuda_stream(handle); - - paramsPCA prms_with_dims = prms; - prms_with_dims.n_rows = static_cast(input.extent(0)); - prms_with_dims.n_cols = static_cast(input.extent(1)); - detail::pcaFitTransform(handle, - input.data_handle(), - trans_input.data_handle(), - components.data_handle(), - explained_var.data_handle(), - explained_var_ratio.data_handle(), - singular_vals.data_handle(), - mu.data_handle(), - noise_vars.data_handle(), - prms_with_dims, - stream, + input, + trans_input, + components, + explained_var, + explained_var_ratio, + singular_vals, + mu, + noise_vars, + prms, flip_signs_based_on_U); } @@ -144,20 +129,7 @@ void pca_inverse_transform(raft::resources const& handle, raft::device_vector_view mu, raft::device_matrix_view output) { - auto stream = resource::get_cuda_stream(handle); - - paramsPCA prms_with_dims = prms; - prms_with_dims.n_rows = static_cast(output.extent(0)); - prms_with_dims.n_cols = static_cast(output.extent(1)); - - detail::pcaInverseTransform(handle, - trans_input.data_handle(), - components.data_handle(), - singular_vals.data_handle(), - mu.data_handle(), - output.data_handle(), - prms_with_dims, - stream); + detail::pcaInverseTransform(handle, trans_input, components, singular_vals, mu, output, prms); } /** @@ -183,20 +155,7 @@ void pca_transform(raft::resources const& handle, raft::device_vector_view mu, raft::device_matrix_view trans_input) { - auto stream = resource::get_cuda_stream(handle); - - paramsPCA prms_with_dims = prms; - prms_with_dims.n_rows = static_cast(input.extent(0)); - prms_with_dims.n_cols = static_cast(input.extent(1)); - - detail::pcaTransform(handle, - input.data_handle(), - components.data_handle(), - trans_input.data_handle(), - singular_vals.data_handle(), - mu.data_handle(), - prms_with_dims, - stream); + detail::pcaTransform(handle, input, components, trans_input, singular_vals, mu, prms); } /** @@ -205,7 +164,7 @@ void pca_transform(raft::resources const& handle, * @tparam math_t data-type upon which the math operation will be performed * @tparam idx_t integer type used for indexing * @param[in] handle raft::resources - * @param[in] prms tSVD parameters (controls n_components, algorithm) + * @param[in] prms tSVD parameters (controls n_components, algorithm). n_rows must be set by caller. * @param[inout] in covariance matrix [n_cols x n_cols] (col-major). Overwritten. * @param[out] components truncated eigenvectors [n_components x n_cols] (col-major) * @param[out] explained_var explained variances [n_components] @@ -221,20 +180,8 @@ void trunc_comp_exp_vars(raft::resources const& handle, raft::device_vector_view explained_var_ratio, raft::device_scalar_view noise_vars) { - auto stream = resource::get_cuda_stream(handle); - - paramsTSVD prms_with_dims = prms; - prms_with_dims.n_rows = static_cast(in.extent(0)); - prms_with_dims.n_cols = static_cast(in.extent(1)); - - detail::truncCompExpVars(handle, - in.data_handle(), - components.data_handle(), - explained_var.data_handle(), - explained_var_ratio.data_handle(), - noise_vars.data_handle(), - prms_with_dims, - stream); + detail::truncCompExpVars( + handle, in, components, explained_var, explained_var_ratio, noise_vars, prms); } /** @} */ // end group pca diff --git a/cpp/include/raft/linalg/tsvd.cuh b/cpp/include/raft/linalg/tsvd.cuh index fe4eb9437b..506539c331 100644 --- a/cpp/include/raft/linalg/tsvd.cuh +++ b/cpp/include/raft/linalg/tsvd.cuh @@ -8,7 +8,6 @@ #include "detail/tsvd.cuh" #include -#include namespace raft::linalg { @@ -37,19 +36,7 @@ void tsvd_fit(raft::resources const& handle, raft::device_vector_view singular_vals, bool flip_signs_based_on_U = false) { - auto stream = resource::get_cuda_stream(handle); - - paramsTSVD prms_with_dims = prms; - prms_with_dims.n_rows = static_cast(input.extent(0)); - prms_with_dims.n_cols = static_cast(input.extent(1)); - - detail::tsvdFit(handle, - input.data_handle(), - components.data_handle(), - singular_vals.data_handle(), - prms_with_dims, - stream, - flip_signs_based_on_U); + detail::tsvdFit(handle, input, components, singular_vals, prms, flip_signs_based_on_U); } /** @@ -81,21 +68,14 @@ void tsvd_fit_transform(raft::resources const& handle, raft::device_vector_view singular_vals, bool flip_signs_based_on_U = false) { - auto stream = resource::get_cuda_stream(handle); - - paramsTSVD prms_with_dims = prms; - prms_with_dims.n_rows = static_cast(input.extent(0)); - prms_with_dims.n_cols = static_cast(input.extent(1)); - detail::tsvdFitTransform(handle, - input.data_handle(), - trans_input.data_handle(), - components.data_handle(), - explained_var.data_handle(), - explained_var_ratio.data_handle(), - singular_vals.data_handle(), - prms_with_dims, - stream, + input, + trans_input, + components, + explained_var, + explained_var_ratio, + singular_vals, + prms, flip_signs_based_on_U); } @@ -118,18 +98,7 @@ void tsvd_transform(raft::resources const& handle, raft::device_matrix_view components, raft::device_matrix_view trans_input) { - auto stream = resource::get_cuda_stream(handle); - - paramsTSVD prms_with_dims = prms; - prms_with_dims.n_rows = static_cast(input.extent(0)); - prms_with_dims.n_cols = static_cast(input.extent(1)); - - detail::tsvdTransform(handle, - input.data_handle(), - components.data_handle(), - trans_input.data_handle(), - prms_with_dims, - stream); + detail::tsvdTransform(handle, input, components, trans_input, prms); } /** @@ -151,18 +120,7 @@ void tsvd_inverse_transform(raft::resources const& handle, raft::device_matrix_view components, raft::device_matrix_view output) { - auto stream = resource::get_cuda_stream(handle); - - paramsTSVD prms_with_dims = prms; - prms_with_dims.n_rows = static_cast(output.extent(0)); - prms_with_dims.n_cols = static_cast(output.extent(1)); - - detail::tsvdInverseTransform(handle, - trans_input.data_handle(), - components.data_handle(), - output.data_handle(), - prms_with_dims, - stream); + detail::tsvdInverseTransform(handle, trans_input, components, output, prms); } /** @@ -183,18 +141,7 @@ void cal_eig(raft::resources const& handle, raft::device_matrix_view components, raft::device_vector_view explained_var) { - auto stream = resource::get_cuda_stream(handle); - - paramsTSVD prms_with_dims = prms; - prms_with_dims.n_rows = static_cast(in.extent(0)); - prms_with_dims.n_cols = static_cast(in.extent(1)); - - detail::calEig(handle, - in.data_handle(), - components.data_handle(), - explained_var.data_handle(), - prms_with_dims, - stream); + detail::calEig(handle, in, components, explained_var, prms); } /** @@ -214,20 +161,7 @@ void sign_flip_components(raft::resources const& handle, bool center, bool flip_signs_based_on_U = false) { - auto stream = resource::get_cuda_stream(handle); - auto n_samples = static_cast(input.extent(0)); - auto n_features = static_cast(input.extent(1)); - auto n_components = static_cast(components.extent(0)); - - detail::signFlipComponents(handle, - input.data_handle(), - components.data_handle(), - n_samples, - n_features, - n_components, - stream, - center, - flip_signs_based_on_U); + detail::signFlipComponents(handle, input, components, center, flip_signs_based_on_U); } /** @} */ // end group tsvd From ac9a88318fdf51b1f6408b916d0c5608a2fa89da Mon Sep 17 00:00:00 2001 From: aamijar Date: Wed, 18 Mar 2026 02:45:24 +0000 Subject: [PATCH 10/12] clean up detail apis 2 --- cpp/include/raft/linalg/detail/pca.cuh | 176 ++++++++++++------------ cpp/include/raft/linalg/detail/tsvd.cuh | 153 ++++++++++---------- cpp/include/raft/linalg/pca.cuh | 50 +++---- cpp/include/raft/linalg/pca_types.hpp | 1 - cpp/include/raft/linalg/tsvd.cuh | 28 ++-- docs/source/cpp_api/linalg_solver.rst | 24 ++++ 6 files changed, 222 insertions(+), 210 deletions(-) diff --git a/cpp/include/raft/linalg/detail/pca.cuh b/cpp/include/raft/linalg/detail/pca.cuh index 6a8ed8c17d..86f18d7dad 100644 --- a/cpp/include/raft/linalg/detail/pca.cuh +++ b/cpp/include/raft/linalg/detail/pca.cuh @@ -28,13 +28,13 @@ namespace raft::linalg::detail { template -void truncCompExpVars(raft::resources const& handle, - raft::device_matrix_view in, - raft::device_matrix_view components, - raft::device_vector_view explained_var, - raft::device_vector_view explained_var_ratio, - raft::device_scalar_view noise_vars, - const paramsTSVD& prms) +void trunc_comp_exp_vars(raft::resources const& handle, + const paramsTSVD& prms, + raft::device_matrix_view in, + raft::device_matrix_view components, + raft::device_vector_view explained_var, + raft::device_vector_view explained_var_ratio, + raft::device_scalar_view noise_vars) { auto stream = resource::get_cuda_stream(handle); @@ -46,13 +46,13 @@ void truncCompExpVars(raft::resources const& handle, rmm::device_uvector explained_var_all(static_cast(n_cols), stream); rmm::device_uvector explained_var_ratio_all(static_cast(n_cols), stream); - detail::calEig( + detail::cal_eig( handle, + prms, in, raft::make_device_matrix_view( components_all.data(), n_cols, n_cols), - raft::make_device_vector_view(explained_var_all.data(), n_cols), - prms); + raft::make_device_vector_view(explained_var_all.data(), n_cols)); raft::matrix::trunc_zero_origin( handle, raft::make_device_matrix_view( @@ -94,9 +94,9 @@ void truncCompExpVars(raft::resources const& handle, } /** - * @brief perform fit operation for the pca. Generates eigenvectors, explained vars, singular vals, - * etc. + * @brief perform fit operation for PCA. * @param[in] handle: raft::resources + * @param[in] prms: PCA parameters (n_components, algorithm, whiten, etc.) * @param[inout] input: the data is fitted to PCA. Size n_rows x n_cols (col-major). * @param[out] components: the principal components. Size n_components x n_cols (col-major). * @param[out] explained_var: explained variances. Size n_components. @@ -104,20 +104,19 @@ void truncCompExpVars(raft::resources const& handle, * @param[out] singular_vals: singular values. Size n_components. * @param[out] mu: mean of all features. Size n_cols. * @param[out] noise_vars: noise variance scalar. - * @param[in] prms: data structure that includes all the parameters from input size to algorithm. * @param[in] flip_signs_based_on_U whether to determine signs by U (true) or V.T (false) */ template -void pcaFit(raft::resources const& handle, - 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, - const paramsPCA& prms, - bool flip_signs_based_on_U = false) +void pca_fit(raft::resources const& handle, + const paramsPCA& prms, + 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) { auto stream = resource::get_cuda_stream(handle); auto cublas_handle = raft::resource::get_cublas_handle(handle); @@ -145,14 +144,14 @@ void pcaFit(raft::resources const& handle, prms_with_rows.n_rows = static_cast(n_rows); prms_with_rows.n_cols = static_cast(n_cols); - detail::truncCompExpVars( + detail::trunc_comp_exp_vars( handle, + prms_with_rows, raft::make_device_matrix_view(cov.data(), n_cols, n_cols), components, explained_var, explained_var_ratio, - noise_vars, - prms_with_rows); + noise_vars); math_t scalar = (n_rows - 1); raft::matrix::weighted_sqrt(handle, @@ -166,28 +165,27 @@ void pcaFit(raft::resources const& handle, raft::stats::meanAdd( input.data_handle(), input.data_handle(), mu.data_handle(), n_cols, n_rows, stream); - detail::signFlipComponents(handle, input, components, true, flip_signs_based_on_U); + detail::sign_flip_components(handle, input, components, true, flip_signs_based_on_U); } /** - * @brief performs transform operation for the pca. Transforms the data to eigenspace. + * @brief performs transform operation for PCA. Transforms the data to eigenspace. * @param[in] handle: raft::resources - * @param[inout] input: the data to transform. Size n_rows x n_cols (col-major). Modified - * temporarily (mean-centered then restored). + * @param[in] prms: PCA parameters (n_components, algorithm, whiten, etc.) + * @param[inout] input: the data to transform. Size n_rows x n_cols (col-major). * @param[in] components: principal components. Size n_components x n_cols (col-major). - * @param[out] trans_input: the transformed data. Size n_rows x n_components (col-major). * @param[in] singular_vals: singular values. Size n_components. * @param[in] mu: mean of features. Size n_cols. - * @param[in] prms: data structure that includes all the parameters from input size to algorithm. + * @param[out] trans_input: the transformed data. Size n_rows x n_components (col-major). */ template -void pcaTransform(raft::resources const& handle, - raft::device_matrix_view input, - raft::device_matrix_view components, - raft::device_matrix_view trans_input, - raft::device_vector_view singular_vals, - raft::device_vector_view mu, - const paramsPCA& prms) +void pca_transform(raft::resources const& handle, + const paramsPCA& prms, + 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 stream = resource::get_cuda_stream(handle); @@ -217,40 +215,39 @@ void pcaTransform(raft::resources const& handle, raft::stats::meanCenter( input.data_handle(), input.data_handle(), mu.data_handle(), n_cols, n_rows, stream); - detail::tsvdTransform(handle, - input, - raft::make_device_matrix_view( - components_copy.data(), n_components, n_cols), - trans_input, - prms); + detail::tsvd_transform(handle, + prms, + input, + raft::make_device_matrix_view( + components_copy.data(), n_components, n_cols), + trans_input); raft::stats::meanAdd( input.data_handle(), input.data_handle(), mu.data_handle(), n_cols, n_rows, stream); } /** - * @brief performs inverse transform operation for the pca. Transforms the transformed data back to - * original data. + * @brief performs inverse transform operation for PCA. * @param[in] handle: raft::resources + * @param[in] prms: PCA parameters (n_components, algorithm, whiten, etc.) * @param[in] trans_input: the transformed data. Size n_rows x n_components (col-major). * @param[in] components: principal components. Size n_components x n_cols (col-major). * @param[in] singular_vals: singular values. Size n_components. * @param[in] mu: mean of features. Size n_cols. - * @param[out] input: the reconstructed data. Size n_rows x n_cols (col-major). - * @param[in] prms: data structure that includes all the parameters from input size to algorithm. + * @param[out] output: the reconstructed data. Size n_rows x n_cols (col-major). */ template -void pcaInverseTransform(raft::resources const& handle, - 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 input, - const paramsPCA& prms) +void pca_inverse_transform(raft::resources const& handle, + const paramsPCA& prms, + 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 stream = resource::get_cuda_stream(handle); - auto n_rows = input.extent(0); - auto n_cols = input.extent(1); + auto n_rows = output.extent(0); + auto n_cols = output.extent(1); auto n_components = components.extent(0); ASSERT(n_cols > 1, "Parameter n_cols: number of columns cannot be less than two"); @@ -274,20 +271,20 @@ void pcaInverseTransform(raft::resources const& handle, n_components)); } - detail::tsvdInverseTransform(handle, - trans_input, - raft::make_device_matrix_view( - components_copy.data(), n_components, n_cols), - input, - prms); + detail::tsvd_inverse_transform(handle, + prms, + trans_input, + raft::make_device_matrix_view( + components_copy.data(), n_components, n_cols), + output); raft::stats::meanAdd( - input.data_handle(), input.data_handle(), mu.data_handle(), n_cols, n_rows, stream); + output.data_handle(), output.data_handle(), mu.data_handle(), n_cols, n_rows, stream); } /** - * @brief perform fit and transform operations for the pca. Generates transformed data, - * eigenvectors, explained vars, singular vals, etc. + * @brief perform fit and transform operations for PCA. * @param[in] handle: raft::resources + * @param[in] prms: PCA parameters (n_components, algorithm, whiten, etc.) * @param[inout] input: the data is fitted to PCA. Size n_rows x n_cols (col-major). * @param[out] trans_input: the transformed data. Size n_rows x n_components (col-major). * @param[out] components: the principal components. Size n_components x n_cols (col-major). @@ -296,33 +293,32 @@ void pcaInverseTransform(raft::resources const& handle, * @param[out] singular_vals: singular values. Size n_components. * @param[out] mu: mean of all features. Size n_cols. * @param[out] noise_vars: noise variance scalar. - * @param[in] prms: data structure that includes all the parameters from input size to algorithm. * @param[in] flip_signs_based_on_U whether to determine signs by U (true) or V.T (false) */ template -void pcaFitTransform(raft::resources const& handle, - 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, - const paramsPCA& prms, - bool flip_signs_based_on_U = false) +void pca_fit_transform(raft::resources const& handle, + const paramsPCA& prms, + 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) { - detail::pcaFit(handle, - input, - components, - explained_var, - explained_var_ratio, - singular_vals, - mu, - noise_vars, - prms, - flip_signs_based_on_U); - detail::pcaTransform(handle, input, components, trans_input, singular_vals, mu, prms); + detail::pca_fit(handle, + prms, + input, + components, + explained_var, + explained_var_ratio, + singular_vals, + mu, + noise_vars, + flip_signs_based_on_U); + detail::pca_transform(handle, prms, input, components, singular_vals, mu, trans_input); } }; // end namespace raft::linalg::detail diff --git a/cpp/include/raft/linalg/detail/tsvd.cuh b/cpp/include/raft/linalg/detail/tsvd.cuh index 5494190e49..1b5d813ac7 100644 --- a/cpp/include/raft/linalg/detail/tsvd.cuh +++ b/cpp/include/raft/linalg/detail/tsvd.cuh @@ -40,13 +40,13 @@ namespace raft::linalg::detail { template -void calCompExpVarsSvd(raft::resources const& handle, - raft::device_matrix_view in, - raft::device_matrix_view components, - raft::device_vector_view singular_vals, - raft::device_vector_view explained_vars, - raft::device_vector_view explained_var_ratio, - const paramsTSVD& prms) +void cal_comp_exp_vars_svd(raft::resources const& handle, + const paramsTSVD& prms, + raft::device_matrix_view in, + raft::device_matrix_view components, + raft::device_vector_view singular_vals, + raft::device_vector_view explained_vars, + raft::device_vector_view explained_var_ratio) { auto stream = resource::get_cuda_stream(handle); auto cusolver_handle = raft::resource::get_cusolver_dn_handle(handle); @@ -103,11 +103,11 @@ void calCompExpVarsSvd(raft::resources const& handle, } template -void calEig(raft::resources const& handle, - raft::device_matrix_view in, - raft::device_matrix_view components, - raft::device_vector_view explained_var, - const paramsTSVD& prms) +void cal_eig(raft::resources const& handle, + const paramsTSVD& prms, + raft::device_matrix_view in, + raft::device_matrix_view components, + raft::device_vector_view explained_var) { auto stream = resource::get_cuda_stream(handle); auto cusolver_handle = raft::resource::get_cusolver_dn_handle(handle); @@ -147,21 +147,19 @@ void calEig(raft::resources const& handle, } /** - * @defgroup sign flip for PCA and tSVD. This is used to stabilize the sign of column major eigen - * vectors + * @brief sign flip for PCA and tSVD. Stabilizes the sign of column major eigenvectors. * @param handle: raft::resources * @param input: input data [n_samples x n_features] (col-major) * @param components: components matrix [n_components x n_features] (col-major) * @param center whether to mean-center input before computing signs * @param flip_signs_based_on_U whether to determine signs by U (true) or V.T (false) - * @{ */ template -void signFlipComponents(raft::resources const& handle, - raft::device_matrix_view input, - raft::device_matrix_view components, - bool center, - bool flip_signs_based_on_U = false) +void sign_flip_components(raft::resources const& handle, + raft::device_matrix_view input, + raft::device_matrix_view components, + bool center, + bool flip_signs_based_on_U = false) { auto stream = resource::get_cuda_stream(handle); auto n_samples = input.extent(0); @@ -240,17 +238,15 @@ void signFlipComponents(raft::resources const& handle, } /** - * @defgroup sign flip for PCA and tSVD. This is used to stabilize the sign of column major eigen - * vectors + * @brief sign flip for PCA and tSVD. Stabilizes the sign of column major eigenvectors. * @param handle: raft::resources * @param input: input matrix [n_rows x n_cols] (col-major). Modified in place. * @param components: components matrix [n_rows x n_cols_comp] (col-major). Modified in place. - * @{ */ template -void signFlip(raft::resources const& handle, - raft::device_matrix_view input, - raft::device_matrix_view components) +void sign_flip(raft::resources const& handle, + raft::device_matrix_view input, + raft::device_matrix_view components) { auto stream = resource::get_cuda_stream(handle); auto n_rows = input.extent(0); @@ -292,22 +288,21 @@ void signFlip(raft::resources const& handle, } /** - * @brief perform fit operation for the tsvd. Generates eigenvectors, explained vars, singular vals, - * etc. + * @brief perform fit operation for the tsvd. * @param[in] handle: raft::resources + * @param[in] prms: data structure that includes all the parameters from input size to algorithm. * @param[in] input: the data is fitted to tSVD. Size n_rows x n_cols (col-major). * @param[out] components: the principal components. Size n_components x n_cols (col-major). * @param[out] singular_vals: singular values of the data. Size n_components. - * @param[in] prms: data structure that includes all the parameters from input size to algorithm. * @param[in] flip_signs_based_on_U whether to determine signs by U (true) or V.T (false) */ template -void tsvdFit(raft::resources const& handle, - raft::device_matrix_view input, - raft::device_matrix_view components, - raft::device_vector_view singular_vals, - const paramsTSVD& prms, - bool flip_signs_based_on_U = false) +void tsvd_fit(raft::resources const& handle, + const paramsTSVD& prms, + raft::device_matrix_view input, + raft::device_matrix_view components, + raft::device_vector_view singular_vals, + bool flip_signs_based_on_U = false) { auto stream = resource::get_cuda_stream(handle); auto cublas_handle = raft::resource::get_cublas_handle(handle); @@ -345,13 +340,13 @@ void tsvdFit(raft::resources const& handle, rmm::device_uvector components_all(len, stream); rmm::device_uvector explained_var_all(static_cast(n_cols), stream); - detail::calEig(handle, - raft::make_device_matrix_view( - input_cross_mult.data(), n_cols, n_cols), - raft::make_device_matrix_view( - components_all.data(), n_cols, n_cols), - raft::make_device_vector_view(explained_var_all.data(), n_cols), - prms); + detail::cal_eig(handle, + prms, + raft::make_device_matrix_view( + input_cross_mult.data(), n_cols, n_cols), + raft::make_device_matrix_view( + components_all.data(), n_cols, n_cols), + raft::make_device_vector_view(explained_var_all.data(), n_cols)); raft::matrix::trunc_zero_origin( handle, @@ -368,28 +363,28 @@ void tsvdFit(raft::resources const& handle, singular_vals.data_handle(), idx_t(1), n_components), raft::make_host_scalar_view(&scalar)); - detail::signFlipComponents(handle, - input, - raft::make_device_matrix_view( - components.data_handle(), n_components, n_cols), - false, - flip_signs_based_on_U); + detail::sign_flip_components(handle, + input, + raft::make_device_matrix_view( + components.data_handle(), n_components, n_cols), + false, + flip_signs_based_on_U); } /** * @brief performs transform operation for the tsvd. Transforms the data to eigenspace. * @param[in] handle raft::resources + * @param[in] prms: data structure that includes all the parameters from input size to algorithm. * @param[in] input: the data to transform. Size n_rows x n_cols (col-major). * @param[in] components: principal components. Size n_components x n_cols (col-major). * @param[out] trans_input: transformed output. Size n_rows x n_components (col-major). - * @param[in] prms: data structure that includes all the parameters from input size to algorithm. */ template -void tsvdTransform(raft::resources const& handle, - raft::device_matrix_view input, - raft::device_matrix_view components, - raft::device_matrix_view trans_input, - const paramsTSVD& prms) +void tsvd_transform(raft::resources const& handle, + const paramsTSVD& prms, + raft::device_matrix_view input, + raft::device_matrix_view components, + raft::device_matrix_view trans_input) { auto stream = resource::get_cuda_stream(handle); @@ -419,25 +414,24 @@ void tsvdTransform(raft::resources const& handle, } /** - * @brief performs inverse transform operation for the tsvd. Transforms the transformed data back to - * original data. + * @brief performs inverse transform operation for the tsvd. * @param[in] handle raft::resources + * @param[in] prms: data structure that includes all the parameters from input size to algorithm. * @param[in] trans_input: the transformed data. Size n_rows x n_components (col-major). * @param[in] components: principal components. Size n_components x n_cols (col-major). - * @param[out] input: reconstructed output. Size n_rows x n_cols (col-major). - * @param[in] prms: data structure that includes all the parameters from input size to algorithm. + * @param[out] output: reconstructed output. Size n_rows x n_cols (col-major). */ template -void tsvdInverseTransform(raft::resources const& handle, - raft::device_matrix_view trans_input, - raft::device_matrix_view components, - raft::device_matrix_view input, - const paramsTSVD& prms) +void tsvd_inverse_transform(raft::resources const& handle, + const paramsTSVD& prms, + raft::device_matrix_view trans_input, + raft::device_matrix_view components, + raft::device_matrix_view output) { auto stream = resource::get_cuda_stream(handle); - auto n_rows = input.extent(0); - auto n_cols = input.extent(1); + auto n_rows = output.extent(0); + auto n_cols = output.extent(1); auto n_components = components.extent(0); ASSERT(n_cols > 1, "Parameter n_cols: number of columns cannot be less than one"); @@ -452,7 +446,7 @@ void tsvdInverseTransform(raft::resources const& handle, n_rows, n_components, components.data_handle(), - input.data_handle(), + output.data_handle(), n_rows, n_cols, CUBLAS_OP_N, @@ -463,28 +457,27 @@ void tsvdInverseTransform(raft::resources const& handle, } /** - * @brief performs fit and transform operations for the tsvd. Generates transformed data, - * eigenvectors, explained vars, singular vals, etc. + * @brief performs fit and transform operations for the tsvd. * @param[in] handle: raft::resources + * @param[in] prms: data structure that includes all the parameters from input size to algorithm. * @param[in] input: the data is fitted to tSVD. Size n_rows x n_cols (col-major). * @param[out] trans_input: the transformed data. Size n_rows x n_components (col-major). * @param[out] components: the principal components. Size n_components x n_cols (col-major). * @param[out] explained_var: explained variances. Size n_components. * @param[out] explained_var_ratio: ratio of explained variance to total. Size n_components. * @param[out] singular_vals: singular values of the data. Size n_components. - * @param[in] prms: data structure that includes all the parameters from input size to algorithm. * @param[in] flip_signs_based_on_U whether to determine signs by U (true) or V.T (false) */ template -void tsvdFitTransform(raft::resources const& handle, - 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, - const paramsTSVD& prms, - bool flip_signs_based_on_U = false) +void tsvd_fit_transform(raft::resources const& handle, + const paramsTSVD& prms, + 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, + bool flip_signs_based_on_U = false) { auto stream = resource::get_cuda_stream(handle); @@ -492,8 +485,8 @@ void tsvdFitTransform(raft::resources const& handle, auto n_cols = input.extent(1); auto n_components = components.extent(0); - detail::tsvdFit(handle, input, components, singular_vals, prms, flip_signs_based_on_U); - detail::tsvdTransform(handle, input, components, trans_input, prms); + detail::tsvd_fit(handle, prms, input, components, singular_vals, flip_signs_based_on_U); + detail::tsvd_transform(handle, prms, input, components, trans_input); rmm::device_uvector mu_trans(static_cast(n_components), stream); raft::stats::mean( diff --git a/cpp/include/raft/linalg/pca.cuh b/cpp/include/raft/linalg/pca.cuh index 5136476884..b9d5a23a5f 100644 --- a/cpp/include/raft/linalg/pca.cuh +++ b/cpp/include/raft/linalg/pca.cuh @@ -47,16 +47,16 @@ void pca_fit(raft::resources const& handle, raft::device_scalar_view noise_vars, bool flip_signs_based_on_U = false) { - detail::pcaFit(handle, - input, - components, - explained_var, - explained_var_ratio, - singular_vals, - mu, - noise_vars, - prms, - flip_signs_based_on_U); + detail::pca_fit(handle, + prms, + input, + components, + explained_var, + explained_var_ratio, + singular_vals, + mu, + noise_vars, + flip_signs_based_on_U); } /** @@ -93,17 +93,17 @@ void pca_fit_transform(raft::resources const& handle, raft::device_scalar_view noise_vars, bool flip_signs_based_on_U = false) { - detail::pcaFitTransform(handle, - input, - trans_input, - components, - explained_var, - explained_var_ratio, - singular_vals, - mu, - noise_vars, - prms, - flip_signs_based_on_U); + detail::pca_fit_transform(handle, + prms, + input, + trans_input, + components, + explained_var, + explained_var_ratio, + singular_vals, + mu, + noise_vars, + flip_signs_based_on_U); } /** @@ -129,7 +129,7 @@ void pca_inverse_transform(raft::resources const& handle, raft::device_vector_view mu, raft::device_matrix_view output) { - detail::pcaInverseTransform(handle, trans_input, components, singular_vals, mu, output, prms); + detail::pca_inverse_transform(handle, prms, trans_input, components, singular_vals, mu, output); } /** @@ -155,7 +155,7 @@ void pca_transform(raft::resources const& handle, raft::device_vector_view mu, raft::device_matrix_view trans_input) { - detail::pcaTransform(handle, input, components, trans_input, singular_vals, mu, prms); + detail::pca_transform(handle, prms, input, components, singular_vals, mu, trans_input); } /** @@ -180,8 +180,8 @@ void trunc_comp_exp_vars(raft::resources const& handle, raft::device_vector_view explained_var_ratio, raft::device_scalar_view noise_vars) { - detail::truncCompExpVars( - handle, in, components, explained_var, explained_var_ratio, noise_vars, prms); + detail::trunc_comp_exp_vars( + handle, prms, in, components, explained_var, explained_var_ratio, noise_vars); } /** @} */ // end group pca diff --git a/cpp/include/raft/linalg/pca_types.hpp b/cpp/include/raft/linalg/pca_types.hpp index 1eba5e7fdb..fe592bb53b 100644 --- a/cpp/include/raft/linalg/pca_types.hpp +++ b/cpp/include/raft/linalg/pca_types.hpp @@ -24,7 +24,6 @@ enum class solver : int { struct paramsTSVD { std::size_t n_rows = 0; std::size_t n_cols = 0; - int gpu_id = 0; float tol = 0.0; uint64_t n_iterations = 15; uint64_t n_components = 1; diff --git a/cpp/include/raft/linalg/tsvd.cuh b/cpp/include/raft/linalg/tsvd.cuh index 506539c331..a71a7f8263 100644 --- a/cpp/include/raft/linalg/tsvd.cuh +++ b/cpp/include/raft/linalg/tsvd.cuh @@ -36,7 +36,7 @@ void tsvd_fit(raft::resources const& handle, raft::device_vector_view singular_vals, bool flip_signs_based_on_U = false) { - detail::tsvdFit(handle, input, components, singular_vals, prms, flip_signs_based_on_U); + detail::tsvd_fit(handle, prms, input, components, singular_vals, flip_signs_based_on_U); } /** @@ -68,15 +68,15 @@ void tsvd_fit_transform(raft::resources const& handle, raft::device_vector_view singular_vals, bool flip_signs_based_on_U = false) { - detail::tsvdFitTransform(handle, - input, - trans_input, - components, - explained_var, - explained_var_ratio, - singular_vals, - prms, - flip_signs_based_on_U); + detail::tsvd_fit_transform(handle, + prms, + input, + trans_input, + components, + explained_var, + explained_var_ratio, + singular_vals, + flip_signs_based_on_U); } /** @@ -98,7 +98,7 @@ void tsvd_transform(raft::resources const& handle, raft::device_matrix_view components, raft::device_matrix_view trans_input) { - detail::tsvdTransform(handle, input, components, trans_input, prms); + detail::tsvd_transform(handle, prms, input, components, trans_input); } /** @@ -120,7 +120,7 @@ void tsvd_inverse_transform(raft::resources const& handle, raft::device_matrix_view components, raft::device_matrix_view output) { - detail::tsvdInverseTransform(handle, trans_input, components, output, prms); + detail::tsvd_inverse_transform(handle, prms, trans_input, components, output); } /** @@ -141,7 +141,7 @@ void cal_eig(raft::resources const& handle, raft::device_matrix_view components, raft::device_vector_view explained_var) { - detail::calEig(handle, in, components, explained_var, prms); + detail::cal_eig(handle, prms, in, components, explained_var); } /** @@ -161,7 +161,7 @@ void sign_flip_components(raft::resources const& handle, bool center, bool flip_signs_based_on_U = false) { - detail::signFlipComponents(handle, input, components, center, flip_signs_based_on_U); + detail::sign_flip_components(handle, input, components, center, flip_signs_based_on_U); } /** @} */ // end group tsvd diff --git a/docs/source/cpp_api/linalg_solver.rst b/docs/source/cpp_api/linalg_solver.rst index 1a811e072a..7d81aa17e0 100644 --- a/docs/source/cpp_api/linalg_solver.rst +++ b/docs/source/cpp_api/linalg_solver.rst @@ -64,3 +64,27 @@ namespace *raft::linalg* :project: RAFT :members: :content-only: + +PCA +--- + +``#include `` + +namespace *raft::linalg* + +.. doxygengroup:: pca + :project: RAFT + :members: + :content-only: + +Truncated SVD +------------- + +``#include `` + +namespace *raft::linalg* + +.. doxygengroup:: tsvd + :project: RAFT + :members: + :content-only: From 202c87c682f34a8765ef123751c49412ac2ce5b9 Mon Sep 17 00:00:00 2001 From: aamijar Date: Wed, 18 Mar 2026 07:37:23 +0000 Subject: [PATCH 11/12] remove n_rows, n_cols, n_components from struct --- cpp/include/raft/linalg/detail/pca.cuh | 23 ++++----- cpp/include/raft/linalg/detail/tsvd.cuh | 9 ++-- cpp/include/raft/linalg/pca.cuh | 8 +-- cpp/include/raft/linalg/pca_types.hpp | 3 -- cpp/tests/linalg/pca.cu | 68 +++++++++++++------------ cpp/tests/linalg/tsvd.cu | 54 ++++++++++---------- 6 files changed, 82 insertions(+), 83 deletions(-) diff --git a/cpp/include/raft/linalg/detail/pca.cuh b/cpp/include/raft/linalg/detail/pca.cuh index 86f18d7dad..d3d34c5549 100644 --- a/cpp/include/raft/linalg/detail/pca.cuh +++ b/cpp/include/raft/linalg/detail/pca.cuh @@ -34,7 +34,8 @@ void trunc_comp_exp_vars(raft::resources const& handle, raft::device_matrix_view components, raft::device_vector_view explained_var, raft::device_vector_view explained_var_ratio, - raft::device_scalar_view noise_vars) + raft::device_scalar_view noise_vars, + std::size_t n_rows) { auto stream = resource::get_cuda_stream(handle); @@ -78,7 +79,7 @@ void trunc_comp_exp_vars(raft::resources const& handle, explained_var_ratio.data_handle(), n_components, idx_t(1))); if (static_cast(n_components) < static_cast(n_cols) && - static_cast(n_components) < prms.n_rows) { + static_cast(n_components) < n_rows) { raft::stats::mean(noise_vars.data_handle(), explained_var_all.data() + static_cast(n_components), std::size_t{1}, @@ -124,13 +125,12 @@ void pca_fit(raft::resources const& handle, auto n_rows = input.extent(0); auto n_cols = input.extent(1); + auto n_components = components.extent(0); + ASSERT(n_cols > 1, "Parameter n_cols: number of columns cannot be less than two"); ASSERT(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 = static_cast(prms.n_components); - if (n_components > n_cols) n_components = n_cols; + ASSERT(n_components > 0, "Parameter n_components: number of components cannot be less than one"); + ASSERT(n_components <= n_cols, "n_components cannot exceed n_cols"); raft::stats::mean(mu.data_handle(), input.data_handle(), n_cols, n_rows, false, stream); @@ -140,18 +140,15 @@ void pca_fit(raft::resources const& handle, raft::stats::cov( handle, cov.data(), input.data_handle(), mu.data_handle(), n_cols, n_rows, true, true, stream); - paramsPCA prms_with_rows = prms; - prms_with_rows.n_rows = static_cast(n_rows); - prms_with_rows.n_cols = static_cast(n_cols); - detail::trunc_comp_exp_vars( handle, - prms_with_rows, + prms, raft::make_device_matrix_view(cov.data(), n_cols, n_cols), components, explained_var, explained_var_ratio, - noise_vars); + noise_vars, + static_cast(n_rows)); math_t scalar = (n_rows - 1); raft::matrix::weighted_sqrt(handle, diff --git a/cpp/include/raft/linalg/detail/tsvd.cuh b/cpp/include/raft/linalg/detail/tsvd.cuh index 1b5d813ac7..8c6ba4cbfe 100644 --- a/cpp/include/raft/linalg/detail/tsvd.cuh +++ b/cpp/include/raft/linalg/detail/tsvd.cuh @@ -310,13 +310,12 @@ void tsvd_fit(raft::resources const& handle, auto n_rows = input.extent(0); auto n_cols = input.extent(1); + auto n_components = components.extent(0); + ASSERT(n_cols > 1, "Parameter n_cols: number of columns cannot be less than two"); ASSERT(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 = static_cast(prms.n_components); - if (n_components > n_cols) n_components = n_cols; + ASSERT(n_components > 0, "Parameter n_components: number of components cannot be less than one"); + ASSERT(n_components <= n_cols, "n_components cannot exceed n_cols"); auto len = static_cast(n_cols * n_cols); rmm::device_uvector input_cross_mult(len, stream); diff --git a/cpp/include/raft/linalg/pca.cuh b/cpp/include/raft/linalg/pca.cuh index b9d5a23a5f..218d12cf77 100644 --- a/cpp/include/raft/linalg/pca.cuh +++ b/cpp/include/raft/linalg/pca.cuh @@ -164,12 +164,13 @@ void pca_transform(raft::resources const& handle, * @tparam math_t data-type upon which the math operation will be performed * @tparam idx_t integer type used for indexing * @param[in] handle raft::resources - * @param[in] prms tSVD parameters (controls n_components, algorithm). n_rows must be set by caller. + * @param[in] prms tSVD parameters (controls algorithm, tolerance, iterations) * @param[inout] in covariance matrix [n_cols x n_cols] (col-major). Overwritten. * @param[out] components truncated eigenvectors [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] noise_vars noise variance scalar + * @param[in] n_rows number of rows in the original data (needed for noise variance computation) */ template void trunc_comp_exp_vars(raft::resources const& handle, @@ -178,10 +179,11 @@ void trunc_comp_exp_vars(raft::resources const& handle, raft::device_matrix_view components, raft::device_vector_view explained_var, raft::device_vector_view explained_var_ratio, - raft::device_scalar_view noise_vars) + raft::device_scalar_view noise_vars, + std::size_t n_rows) { detail::trunc_comp_exp_vars( - handle, prms, in, components, explained_var, explained_var_ratio, noise_vars); + handle, prms, in, components, explained_var, explained_var_ratio, noise_vars, n_rows); } /** @} */ // end group pca diff --git a/cpp/include/raft/linalg/pca_types.hpp b/cpp/include/raft/linalg/pca_types.hpp index fe592bb53b..cca3b09155 100644 --- a/cpp/include/raft/linalg/pca_types.hpp +++ b/cpp/include/raft/linalg/pca_types.hpp @@ -22,11 +22,8 @@ enum class solver : int { /** @brief Parameters for TSVD (and base for PCA). */ struct paramsTSVD { - std::size_t n_rows = 0; - std::size_t n_cols = 0; float tol = 0.0; uint64_t n_iterations = 15; - uint64_t n_components = 1; solver algorithm = solver::COV_EIG_DQ; }; diff --git a/cpp/tests/linalg/pca.cu b/cpp/tests/linalg/pca.cu index beb6220432..a1d80953e2 100644 --- a/cpp/tests/linalg/pca.cu +++ b/cpp/tests/linalg/pca.cu @@ -73,7 +73,8 @@ class PcaTest : public ::testing::TestWithParam> { 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; + int n_components = params.n_col; + int len_comp = params.n_col * params.n_col; rmm::device_uvector explained_var_ratio(params.n_col, stream); rmm::device_uvector singular_vals(params.n_col, stream); rmm::device_uvector mean(params.n_col, stream); @@ -88,27 +89,27 @@ class PcaTest : public ::testing::TestWithParam> { raft::update_device( explained_vars_ref.data(), explained_vars_ref_h.data(), params.n_col, stream); + std::size_t n_rows = params.n_row; + std::size_t n_cols = params.n_col; + paramsPCA prms; - prms.n_cols = params.n_col; - prms.n_rows = params.n_row; - prms.n_components = params.n_col; - prms.whiten = false; + prms.whiten = false; if (params.algo == 0) prms.algorithm = solver::COV_EIG_DQ; else prms.algorithm = solver::COV_EIG_JACOBI; - auto input_view = raft::make_device_matrix_view( - data.data(), prms.n_rows, prms.n_cols); + auto input_view = + raft::make_device_matrix_view(data.data(), n_rows, n_cols); auto components_view = raft::make_device_matrix_view( - components.data(), prms.n_components, prms.n_cols); + components.data(), n_components, n_cols); auto explained_var_view = - raft::make_device_vector_view(explained_vars.data(), prms.n_components); + raft::make_device_vector_view(explained_vars.data(), n_components); auto explained_var_ratio_view = - raft::make_device_vector_view(explained_var_ratio.data(), prms.n_components); + raft::make_device_vector_view(explained_var_ratio.data(), n_components); auto singular_vals_view = - raft::make_device_vector_view(singular_vals.data(), prms.n_components); - auto mu_view = raft::make_device_vector_view(mean.data(), prms.n_cols); + raft::make_device_vector_view(singular_vals.data(), n_components); + auto mu_view = raft::make_device_vector_view(mean.data(), n_cols); auto noise_vars_view = raft::make_device_scalar_view(noise_vars.data()); pca_fit(handle, @@ -122,13 +123,13 @@ class PcaTest : public ::testing::TestWithParam> { noise_vars_view); auto trans_data_view = raft::make_device_matrix_view( - trans_data.data(), prms.n_rows, prms.n_components); + trans_data.data(), n_rows, n_components); pca_transform( handle, prms, input_view, components_view, singular_vals_view, mu_view, trans_data_view); auto data_back_view = raft::make_device_matrix_view( - data_back.data(), prms.n_rows, prms.n_cols); + data_back.data(), n_rows, n_cols); pca_inverse_transform( handle, prms, trans_data_view, components_view, singular_vals_view, mu_view, data_back_view); @@ -139,40 +140,41 @@ class PcaTest : public ::testing::TestWithParam> { raft::random::Rng r(params.seed, raft::random::GenPC); int len = params.len2; + std::size_t n_rows = params.n_row2; + std::size_t n_cols = params.n_col2; + std::size_t n_components = params.n_col2; + paramsPCA prms; - prms.n_cols = params.n_col2; - prms.n_rows = params.n_row2; - prms.n_components = params.n_col2; - prms.whiten = false; + prms.whiten = false; if (params.algo == 0) prms.algorithm = solver::COV_EIG_DQ; else if (params.algo == 1) prms.algorithm = solver::COV_EIG_JACOBI; r.uniform(data2.data(), len, T(-1.0), T(1.0), stream); - rmm::device_uvector data2_trans(prms.n_rows * prms.n_components, stream); + rmm::device_uvector data2_trans(n_rows * n_components, stream); - int len_comp = params.n_col2 * prms.n_components; + int len_comp = params.n_col2 * n_components; rmm::device_uvector components2(len_comp, stream); - rmm::device_uvector explained_vars2(prms.n_components, stream); - rmm::device_uvector explained_var_ratio2(prms.n_components, stream); - rmm::device_uvector singular_vals2(prms.n_components, stream); - rmm::device_uvector mean2(prms.n_cols, stream); + rmm::device_uvector explained_vars2(n_components, stream); + rmm::device_uvector explained_var_ratio2(n_components, stream); + rmm::device_uvector singular_vals2(n_components, stream); + rmm::device_uvector mean2(n_cols, stream); rmm::device_uvector noise_vars2(1, stream); - auto input_view = raft::make_device_matrix_view( - data2.data(), prms.n_rows, prms.n_cols); + auto input_view = + raft::make_device_matrix_view(data2.data(), n_rows, n_cols); auto trans_view = raft::make_device_matrix_view( - data2_trans.data(), prms.n_rows, prms.n_components); + data2_trans.data(), n_rows, n_components); auto comp_view = raft::make_device_matrix_view( - components2.data(), prms.n_components, prms.n_cols); + components2.data(), n_components, n_cols); auto ev_view = - raft::make_device_vector_view(explained_vars2.data(), prms.n_components); + raft::make_device_vector_view(explained_vars2.data(), n_components); auto evr_view = - raft::make_device_vector_view(explained_var_ratio2.data(), prms.n_components); + raft::make_device_vector_view(explained_var_ratio2.data(), n_components); auto sv_view = - raft::make_device_vector_view(singular_vals2.data(), prms.n_components); - auto mu_view = raft::make_device_vector_view(mean2.data(), prms.n_cols); + raft::make_device_vector_view(singular_vals2.data(), n_components); + auto mu_view = raft::make_device_vector_view(mean2.data(), n_cols); auto noise_view = raft::make_device_scalar_view(noise_vars2.data()); pca_fit_transform(handle, @@ -187,7 +189,7 @@ class PcaTest : public ::testing::TestWithParam> { noise_view); auto data2_back_view = raft::make_device_matrix_view( - data2_back.data(), prms.n_rows, prms.n_cols); + data2_back.data(), n_rows, n_cols); pca_inverse_transform(handle, prms, trans_view, comp_view, sv_view, mu_view, data2_back_view); } diff --git a/cpp/tests/linalg/tsvd.cu b/cpp/tests/linalg/tsvd.cu index bfc211df66..9fca1423bf 100644 --- a/cpp/tests/linalg/tsvd.cu +++ b/cpp/tests/linalg/tsvd.cu @@ -74,21 +74,22 @@ class TsvdTest : public ::testing::TestWithParam> { components_ref.resize(len_comp, stream); raft::update_device(components_ref.data(), components_ref_h.data(), len_comp, stream); + std::size_t n_rows = params.n_row; + std::size_t n_cols = params.n_col; + std::size_t n_components = params.n_col; + paramsTSVD prms; - prms.n_cols = params.n_col; - prms.n_rows = params.n_row; - prms.n_components = params.n_col; if (params.algo == 0) prms.algorithm = solver::COV_EIG_DQ; else prms.algorithm = solver::COV_EIG_JACOBI; - auto input_view = raft::make_device_matrix_view( - data.data(), prms.n_rows, prms.n_cols); + auto input_view = + raft::make_device_matrix_view(data.data(), n_rows, n_cols); auto components_view = raft::make_device_matrix_view( - components.data(), prms.n_components, prms.n_cols); + components.data(), n_components, n_cols); auto singular_vals_view = - raft::make_device_vector_view(singular_vals.data(), prms.n_components); + raft::make_device_vector_view(singular_vals.data(), n_components); tsvd_fit(handle, prms, input_view, components_view, singular_vals_view); } @@ -98,16 +99,17 @@ class TsvdTest : public ::testing::TestWithParam> { raft::random::Rng r(params.seed, raft::random::GenPC); int len = params.n_row2 * params.n_col2; + std::size_t n_rows = params.n_row2; + std::size_t n_cols = params.n_col2; + std::size_t n_components = params.n_col2; + paramsTSVD prms; - prms.n_cols = params.n_col2; - prms.n_rows = params.n_row2; - prms.n_components = params.n_col2; if (params.algo == 0) prms.algorithm = solver::COV_EIG_DQ; else if (params.algo == 1) prms.algorithm = solver::COV_EIG_JACOBI; else - prms.n_components = params.n_col2 - 15; + n_components = params.n_col2 - 15; data2.resize(len, stream); int redundant_cols = int(params.redundancy * params.n_col2); @@ -122,37 +124,37 @@ class TsvdTest : public ::testing::TestWithParam> { redundant_len * sizeof(T), cudaMemcpyDeviceToDevice, stream)); - rmm::device_uvector data2_trans(prms.n_rows * prms.n_components, stream); + rmm::device_uvector data2_trans(n_rows * n_components, stream); - int len_comp = params.n_col2 * prms.n_components; + int len_comp = params.n_col2 * n_components; rmm::device_uvector components2(len_comp, stream); - rmm::device_uvector explained_vars2(prms.n_components, stream); - rmm::device_uvector explained_var_ratio2(prms.n_components, stream); - rmm::device_uvector singular_vals2(prms.n_components, stream); + rmm::device_uvector explained_vars2(n_components, stream); + rmm::device_uvector explained_var_ratio2(n_components, stream); + rmm::device_uvector singular_vals2(n_components, stream); - auto input_view = raft::make_device_matrix_view( - data2.data(), prms.n_rows, prms.n_cols); + auto input_view = + raft::make_device_matrix_view(data2.data(), n_rows, n_cols); auto trans_view = raft::make_device_matrix_view( - data2_trans.data(), prms.n_rows, prms.n_components); + data2_trans.data(), n_rows, n_components); auto comp_view = raft::make_device_matrix_view( - components2.data(), prms.n_components, prms.n_cols); + components2.data(), n_components, n_cols); auto ev_view = - raft::make_device_vector_view(explained_vars2.data(), prms.n_components); + raft::make_device_vector_view(explained_vars2.data(), n_components); auto evr_view = - raft::make_device_vector_view(explained_var_ratio2.data(), prms.n_components); + raft::make_device_vector_view(explained_var_ratio2.data(), n_components); auto sv_view = - raft::make_device_vector_view(singular_vals2.data(), prms.n_components); + raft::make_device_vector_view(singular_vals2.data(), n_components); tsvd_fit_transform(handle, prms, input_view, trans_view, comp_view, ev_view, evr_view, sv_view); data2_back.resize(len, stream); auto trans_in_view = raft::make_device_matrix_view( - data2_trans.data(), prms.n_rows, prms.n_components); + data2_trans.data(), n_rows, n_components); auto comp_in_view = raft::make_device_matrix_view( - components2.data(), prms.n_components, prms.n_cols); + components2.data(), n_components, n_cols); auto output_view = raft::make_device_matrix_view( - data2_back.data(), prms.n_rows, prms.n_cols); + data2_back.data(), n_rows, n_cols); tsvd_inverse_transform(handle, prms, trans_in_view, comp_in_view, output_view); } From d97f67d3332d1e25f28eb11cfadbc9c4705ce028 Mon Sep 17 00:00:00 2001 From: Anupam <54245698+aamijar@users.noreply.github.com> Date: Wed, 18 Mar 2026 14:49:31 -0700 Subject: [PATCH 12/12] Update cpp/include/raft/linalg/tsvd.cuh Co-authored-by: Jinsol Park --- cpp/include/raft/linalg/tsvd.cuh | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/cpp/include/raft/linalg/tsvd.cuh b/cpp/include/raft/linalg/tsvd.cuh index a71a7f8263..8f946d5b39 100644 --- a/cpp/include/raft/linalg/tsvd.cuh +++ b/cpp/include/raft/linalg/tsvd.cuh @@ -109,7 +109,7 @@ void tsvd_transform(raft::resources const& handle, * @param[in] handle raft::resources * @param[in] prms data structure that includes all the parameters from input size to algorithm. * @param[in] trans_input the transformed data. Size n_rows x n_components (col-major). - * @param[in] components transpose of the principal components. Size n_components x n_cols + * @param[in] components the principal components. Size n_components x n_cols * (col-major). * @param[out] output the reconstructed data. Size n_rows x n_cols (col-major). */