Update oneDAL PCA infer method

cpp/oneapi/dal/algo/pca/backend/cpu/infer_kernel.cpp

@@ -51,18 +51,20 @@ static result_t call_daal_kernel(const context_cpu& ctx,
 
  const auto daal_data = interop::convert_to_daal_table<Float>(data);
  const auto daal_eigenvectors = interop::convert_to_daal_table<Float>(model.get_eigenvectors());
+ const auto daal_means = interop::convert_to_daal_table<Float>(model.get_means());
+ const auto daal_eigenvalues = interop::convert_to_daal_table<Float>(model.get_eigenvalues());
+
  const auto daal_result =
  interop::convert_to_daal_homogen_table(arr_result, row_count, component_count);
 
  interop::status_to_exception(
  interop::call_daal_kernel<Float, daal_pca_transform_kernel_t>(ctx,
  *daal_data,
  *daal_eigenvectors,
+ daal_means.get(),
  nullptr,
- nullptr,
- nullptr,
+ daal_eigenvalues.get(),
  *daal_result));
-
  return result_t{}.set_transformed_data(
  dal::detail::homogen_table_builder{}.reset(arr_result, row_count, component_count).build());
 }

cpp/oneapi/dal/algo/pca/backend/gpu/infer_kernel_dpc.cpp

@@ -39,25 +39,75 @@ static result_t infer(const context_gpu& ctx, const descriptor_t& desc, const in
  const auto data = input.get_data();
  auto model = input.get_model();
  auto eigenvectors = model.get_eigenvectors();
+ auto means = model.get_means();
+ auto eigenvalues = model.get_eigenvalues();
+
  const std::int64_t row_count = data.get_row_count();
+ ONEDAL_ASSERT(row_count > 0);
+ const std::int64_t col_count = data.get_column_count();
+ ONEDAL_ASSERT(col_count > 0);
  const std::int64_t component_count = get_component_count(desc, data);
+ ONEDAL_ASSERT(component_count > 0);
  dal::detail::check_mul_overflow(row_count, component_count);
 
  const auto data_nd = pr::table2ndarray<Float>(queue, data, sycl::usm::alloc::device);
+ auto means_nd = pr::table2ndarray_1d<Float>(queue, means, sycl::usm::alloc::device);
+ auto mean_centered_data_nd =
+ pr::ndarray<Float, 2>::empty(queue, { row_count, col_count }, sycl::usm::alloc::device);
  const auto eigenvectors_nd =
  pr::table2ndarray<Float>(queue, eigenvectors, sycl::usm::alloc::device);
 
+ sycl::event mean_center_event;
+ {
+ ONEDAL_PROFILER_TASK(elementwise_difference, queue);
+ mean_center_event =
+ pr::elementwise_difference(queue, row_count, data_nd, means_nd, mean_centered_data_nd);
+ }
+
  auto res_nd = pr::ndarray<Float, 2>::empty(queue,
  { row_count, component_count },
  sycl::usm::alloc::device);
  sycl::event gemm_event;
  {
  ONEDAL_PROFILER_TASK(gemm, queue);
- gemm_event = pr::gemm(queue, data_nd, eigenvectors_nd.t(), res_nd, Float(1.0), Float(0.0));
+ gemm_event = pr::gemm(queue,
+ mean_centered_data_nd,
+ eigenvectors_nd.t(),
+ res_nd,
+ Float(1.0),
+ Float(0.0),
+ { mean_center_event });
+ gemm_event.wait_and_throw();
  }
- const auto res_array = res_nd.flatten(queue, { gemm_event });
- auto res_table = homogen_table::wrap(res_array, row_count, component_count);
 
+ if (eigenvalues.has_data()) {
+ auto eigenvalues_nd =
+ pr::table2ndarray_1d<Float>(queue, eigenvalues, sycl::usm::alloc::device);
+ auto sqrt_eigenvalues_nd =
+ pr::ndarray<Float, 1>::empty(queue, { component_count }, sycl::usm::alloc::device);
+ pr::elementwise_sqrt(queue, eigenvalues_nd, sqrt_eigenvalues_nd);
+
+ auto result_whitened_nd = pr::ndarray<Float, 2>::empty(queue,
+ { row_count, component_count },
+ sycl::usm::alloc::device);
+ sycl::event whiten_event;
+ {
+ ONEDAL_PROFILER_TASK(elementwise_division, queue);
+ whiten_event = pr::elementwise_division(queue,
+ row_count,
+ res_nd,
+ sqrt_eigenvalues_nd,
+ result_whitened_nd);
+ }
+
+ const auto res_array_whitened = result_whitened_nd.flatten(queue, { whiten_event });
+ const auto res_table = homogen_table::wrap(res_array_whitened, row_count, component_count);
+
+ return result_t{}.set_transformed_data(res_table);
+ }
+
+ const auto res_array = res_nd.flatten(queue, { gemm_event });
+ const auto res_table = homogen_table::wrap(res_array, row_count, component_count);
  return result_t{}.set_transformed_data(res_table);
 }
 

cpp/oneapi/dal/algo/pca/common.cpp

@@ -61,13 +61,22 @@ template <typename Task>
 class model_impl : public ONEDAL_SERIALIZABLE(pca_dim_reduction_model_impl_id) {
 public:
  table eigenvectors;
+ table pMeans;
+ table pVariances;
+ table eigenvalues;
 
  void serialize(dal::detail::output_archive& ar) const override {
  ar(eigenvectors);
+ ar(pMeans);
+ ar(pVariances);
+ ar(eigenvalues);
  }
 
  void deserialize(dal::detail::input_archive& ar) override {
  ar(eigenvectors);
+ ar(pMeans);
+ ar(pVariances);
+ ar(eigenvalues);
  }
 };
 
@@ -133,6 +142,36 @@ void model<Task>::set_eigenvectors_impl(const table& value) {
  impl_->eigenvectors = value;
 }
 
+template <typename Task>
+const table& model<Task>::get_means() const {
+ return impl_->pMeans;
+}
+
+template <typename Task>
+void model<Task>::set_means_impl(const table& value) {
+ impl_->pMeans = value;
+}
+
+template <typename Task>
+const table& model<Task>::get_variances() const {
+ return impl_->pVariances;
+}
+
+template <typename Task>
+void model<Task>::set_variances_impl(const table& value) {
+ impl_->pVariances = value;
+}
+
+template <typename Task>
+const table& model<Task>::get_eigenvalues() const {
+ return impl_->eigenvalues;
+}
+
+template <typename Task>
+void model<Task>::set_eigenvalues_impl(const table& value) {
+ impl_->eigenvalues = value;
+}
+
 template <typename Task>
 void model<Task>::serialize(dal::detail::output_archive& ar) const {
  dal::detail::serialize_polymorphic_shared(impl_, ar);

cpp/oneapi/dal/algo/pca/common.hpp

@@ -243,8 +243,32 @@ class model : public base {
  return *this;
  }
 
+ const table& get_means() const;
+
+ auto& set_means(const table& value) {
+ set_means_impl(value);
+ return *this;
+ }
+
+ const table& get_variances() const;
+
+ auto& set_variances(const table& value) {
+ set_variances_impl(value);
+ return *this;
+ }
+
+ const table& get_eigenvalues() const;
+
+ auto& set_eigenvalues(const table& value) {
+ set_eigenvalues_impl(value);
+ return *this;
+ }
+
 protected:
  void set_eigenvectors_impl(const table&);
+ void set_means_impl(const table&);
+ void set_variances_impl(const table&);
+ void set_eigenvalues_impl(const table&);
 
 private:
  void serialize(dal::detail::output_archive& ar) const;

cpp/oneapi/dal/backend/primitives/stat/cov.hpp

@@ -37,6 +37,53 @@ sycl::event means(sycl::queue& queue,
  ndview<Float, 1>& means,
  const event_vector& deps = {});
 
+/// Subtract 1-d array from 2-d array
+///
+/// @tparam Float Floating-point type used to perform computations
+///
+/// @param[in] queue The queue
+/// @param[in] row_count The number of rows
+/// @param[in] minuend A 2-d array, from which another array will be subtracted
+/// @param[in] subtrahend A 1-d array, which is to be subtracted from minuend
+/// @param[out] difference The difference between minuend and
+template <typename Float>
+sycl::event elementwise_difference(sycl::queue& queue,
+ std::int64_t row_count,
+ const ndview<Float, 2>& minuend,
+ const ndview<Float, 1>& subtrahend,
+ ndview<Float, 2>& difference,
+ const event_vector& deps = {});
+
+/// SQRT of 1-d array
+///
+/// @tparam Float Floating-point type used to perform computations
+///
+/// @param[in] queue The queue
+/// @param[in] src source array
+/// @param[out] dst destination array with elementwise square-root of source array
+template <typename Float>
+sycl::event elementwise_sqrt(sycl::queue& queue,
+ const ndview<Float, 1>& src,
+ ndview<Float, 1>& dst,
+ const event_vector& deps = {});
+
+/// Divide 2-d array by 1-d array
+///
+/// @tparam Float Floating-point type used to perform computations
+///
+/// @param[in] queue The queue
+/// @param[in] row_count The number of rows
+/// @param[in] numerator A 2-d array
+/// @param[in] denominator A 1-d array
+/// @param[out] quotient The result of the division
+template <typename Float>
+sycl::event elementwise_division(sycl::queue& queue,
+ std::int64_t row_count,
+ const ndview<Float, 2>& numerator,
+ const ndview<Float, 1>& denominator,
+ ndview<Float, 2>& quotient,
+ const event_vector& deps = {});
+
 /// Computes covariance matrix
 ///
 /// @tparam Float Floating-point type used to perform computations

cpp/oneapi/dal/backend/primitives/stat/cov_dpc.cpp

@@ -51,6 +51,98 @@ sycl::event means(sycl::queue& q,
  });
 }
 
+// Subtract 1d array from 2d array elementwise
+template <typename Float>
+sycl::event elementwise_difference(sycl::queue& q,
+ std::int64_t row_count,
+ const ndview<Float, 2>& minuend,
+ const ndview<Float, 1>& subtrahend,
+ ndview<Float, 2>& difference,
+ const event_vector& deps) {
+ ONEDAL_ASSERT(minuend.has_data());
+ ONEDAL_ASSERT(subtrahend.has_data());
+ ONEDAL_ASSERT(difference.has_mutable_data());
+ ONEDAL_ASSERT(is_known_usm(q, minuend.get_data()));
+ ONEDAL_ASSERT(is_known_usm(q, subtrahend.get_data()));
+ ONEDAL_ASSERT(is_known_usm(q, difference.get_mutable_data()));
+ ONEDAL_ASSERT(minuend.get_dimension(0) == difference.get_dimension(0));
+ ONEDAL_ASSERT(minuend.get_dimension(1) == subtrahend.get_dimension(0));
+
+ const auto column_count = minuend.get_dimension(1);
+
+ const Float* minuend_ptr = minuend.get_data();
+ const Float* subtrahend_ptr = subtrahend.get_data();
+ Float* difference_ptr = difference.get_mutable_data();
+
+ return q.submit([&](sycl::handler& cgh) {
+ const auto range = make_range_2d(row_count, column_count);
+ cgh.depends_on(deps);
+ cgh.parallel_for(range, [=](sycl::id<2> id) {
+ const std::int64_t i = id[0];
+ const std::int64_t j = id[1];
+ difference_ptr[i * column_count + j] =
+ minuend_ptr[i * column_count + j] - subtrahend_ptr[j];
+ });
+ });
+}
+
+// SQRT of 1d array elementwise
+template <typename Float>
+sycl::event elementwise_sqrt(sycl::queue& q,
+ const ndview<Float, 1>& src,
+ ndview<Float, 1>& dst,
+ const event_vector& deps) {
+ ONEDAL_ASSERT(src.has_data());
+ ONEDAL_ASSERT(dst.has_mutable_data());
+ ONEDAL_ASSERT(dst.get_shape() == src.get_shape());
+ Float* dst_ptr = dst.get_mutable_data();
+ const Float* src_ptr = src.get_data();
+ const auto num_elems = src.get_dimension(0);
+ return q.submit([&](sycl::handler& cgh) {
+ const auto range = make_range_1d(num_elems);
+ cgh.depends_on(deps);
+ cgh.parallel_for(range, [=](sycl::id<1> id) {
+ dst_ptr[id[0]] = sycl::sqrt(src_ptr[id[0]]);
+ });
+ });
+}
+
+// Divide 2d array by 1d array elementwise
+template <typename Float>
+sycl::event elementwise_division(sycl::queue& q,
+ std::int64_t row_count,
+ const ndview<Float, 2>& numerator,
+ const ndview<Float, 1>& denominator,
+ ndview<Float, 2>& quotient,
+ const event_vector& deps) {
+ ONEDAL_ASSERT(numerator.has_data());
+ ONEDAL_ASSERT(denominator.has_data());
+ ONEDAL_ASSERT(quotient.has_mutable_data());
+ ONEDAL_ASSERT(is_known_usm(q, numerator.get_data()));
+ ONEDAL_ASSERT(is_known_usm(q, denominator.get_data()));
+ ONEDAL_ASSERT(is_known_usm(q, quotient.get_mutable_data()));
+ ONEDAL_ASSERT(numerator.get_dimension(0) == quotient.get_dimension(0));
+ ONEDAL_ASSERT(numerator.get_dimension(1) == denominator.get_dimension(0));
+
+ const auto column_count = numerator.get_dimension(1);
+
+ const Float* numerator_ptr = numerator.get_data();
+ const Float* denominator_ptr = denominator.get_data();
+ Float* quotient_ptr = quotient.get_mutable_data();
+
+ return q.submit([&](sycl::handler& cgh) {
+ const auto range = make_range_2d(row_count, column_count);
+ cgh.depends_on(deps);
+ cgh.parallel_for(range, [=](sycl::id<2> id) {
+ const std::int64_t i = id[0];
+ const std::int64_t j = id[1];
+ ONEDAL_ASSERT(Float(denominator_ptr[j]) != Float(0.0));
+ quotient_ptr[i * column_count + j] =
+ numerator_ptr[i * column_count + j] / denominator_ptr[j];
+ });
+ });
+}
+
 template <typename Float>
 inline sycl::event compute_covariance(sycl::queue& q,
  std::int64_t row_count,
@@ -128,6 +220,7 @@ sycl::event variances(sycl::queue& q,
  });
  });
 }
+
 template <typename Float>
 inline sycl::event prepare_correlation(sycl::queue& q,
  std::int64_t row_count,
@@ -342,6 +435,37 @@ sycl::event correlation_from_covariance(sycl::queue& q,
 INSTANTIATE_MEANS(float)
 INSTANTIATE_MEANS(double)
 
+#define INSTANTIATE_ELEMENWISE_DIFFERENCE(F) \
+ template ONEDAL_EXPORT sycl::event elementwise_difference<F>(sycl::queue&, \
+ std::int64_t, \
+ const ndview<F, 2>&, \
+ const ndview<F, 1>&, \
+ ndview<F, 2>&, \
+ const event_vector&);
+
+INSTANTIATE_ELEMENWISE_DIFFERENCE(float)
+INSTANTIATE_ELEMENWISE_DIFFERENCE(double)
+
+#define INSTANTIATE_ELEMENWISE_DIVISION(F) \
+ template ONEDAL_EXPORT sycl::event elementwise_division<F>(sycl::queue&, \
+ std::int64_t, \
+ const ndview<F, 2>&, \
+ const ndview<F, 1>&, \
+ ndview<F, 2>&, \
+ const event_vector&);
+
+INSTANTIATE_ELEMENWISE_DIVISION(float)
+INSTANTIATE_ELEMENWISE_DIVISION(double)
+
+#define INSTANTIATE_ELEMENWISE_SQRT(F) \
+ template ONEDAL_EXPORT sycl::event elementwise_sqrt<F>(sycl::queue&, \
+ const ndview<F, 1>&, \
+ ndview<F, 1>&, \
+ const event_vector&);
+
+INSTANTIATE_ELEMENWISE_SQRT(float)
+INSTANTIATE_ELEMENWISE_SQRT(double)
+
 #define INSTANTIATE_COV(F) \
  template ONEDAL_EXPORT sycl::event covariance<F>(sycl::queue&, \
  std::int64_t, \