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,        \