feature: add noise_variance computation on oneDAL side

cpp/daal/src/algorithms/pca/pca_dense_base.h

@@ -44,7 +44,9 @@ class PCADenseBase : public Kernel
     services::Status copyTable(NumericTable & source, NumericTable & dest) const;
     services::Status computeExplainedVariancesRatio(const data_management::NumericTable & eigenvalues,
                                                     const data_management::NumericTable & variances,
-                                                    data_management::NumericTable & explained_variances_ratio);
+                                                    data_management::NumericTable & explainedVariancesRatio);
+    services::Status computeNoiseVariances(const data_management::NumericTable & eigenvalues, const data_management::NumericTable & variances,
+                                           double & noiseVariance);
 
 private:
     void signFlipArray(size_t size, algorithmFPType * source) const;

cpp/daal/src/algorithms/pca/pca_dense_base_impl.i

@@ -60,6 +60,35 @@ services::Status PCADenseBase<algorithmFPType, cpu>::computeExplainedVariancesRa
     return services::Status();
 }
 
+template <typename algorithmFPType, CpuType cpu>
+services::Status PCADenseBase<algorithmFPType, cpu>::computeNoiseVariances(const data_management::NumericTable & eigenvalues,
+                                                                           const data_management::NumericTable & variances, double & noiseVariance)
+{
+    const size_t nComponents = eigenvalues.getNumberOfColumns();
+    const size_t nColumns    = variances.getNumberOfColumns();
+
+    ReadRows<algorithmFPType, cpu> eigenValuesBlock(const_cast<data_management::NumericTable &>(eigenvalues), 0, 1);
+    DAAL_CHECK_BLOCK_STATUS(eigenValuesBlock);
+    const algorithmFPType * const eigenValuesArray = eigenValuesBlock.get();
+    ReadRows<algorithmFPType, cpu> variancesBlock(const_cast<data_management::NumericTable &>(variances), 0, 1);
+    DAAL_CHECK_BLOCK_STATUS(variancesBlock);
+    const algorithmFPType * const variancesBlockArray = variancesBlock.get();
+    double totalVariance                              = 0.0;
+    for (size_t i = 0; i < nColumns; i++)
+    {
+        totalVariance += variancesBlockArray[i];
+    }
+
+    double explainedVariance = 0.0;
+    for (size_t i = 0; i < nComponents; i++)
+    {
+        explainedVariance += eigenValuesArray[i];
+    }
+
+    noiseVariance = (totalVariance - explainedVariance) / (nColumns - nComponents);
+    return services::Status();
+}
+
 template <typename algorithmFPType, CpuType cpu>
 services::Status PCADenseBase<algorithmFPType, cpu>::copyTable(NumericTable & source, NumericTable & dest) const
 {

cpp/daal/src/algorithms/pca/pca_dense_correlation_base.h

@@ -46,7 +46,10 @@ class PCACorrelationBase : public PCACorrelationBaseIface<algorithmFPType>, publ
     services::Status computeCorrelationEigenvalues(const data_management::NumericTable & correlation, data_management::NumericTable & eigenvectors,
                                                    data_management::NumericTable & eigenvalues) DAAL_C11_OVERRIDE;
     services::Status computeEigenvectorsInplace(size_t nFeatures, algorithmFPType * eigenvectors, algorithmFPType * eigenvalues);
-    services::Status sortEigenvectorsDescending(size_t nFeatures, algorithmFPType * eigenvectors, algorithmFPType * eigenvalues);
+    services::Status computeEigenvectorsInplaceSyevr(size_t nFeatures, size_t nComponents, algorithmFPType * eigenvectors,
+                                                     algorithmFPType * eigenvalues);
+    DAAL_DEPRECATED services::Status sortEigenvectorsDescending(size_t nFeatures, algorithmFPType * eigenvectors, algorithmFPType * eigenvalues);
+    services::Status sortEigenvectorsDescending(size_t nFeatures, size_t nComponents, algorithmFPType * eigenvectors, algorithmFPType * eigenvalues);
     services::Status computeSingularValues(const data_management::NumericTable & eigenvalues, data_management::NumericTable & variances,
                                            size_t nRows);
     services::Status computeVariancesFromCov(const data_management::NumericTable & correlation, data_management::NumericTable & variances);

cpp/daal/src/algorithms/pca/pca_dense_correlation_base_impl.i

@@ -171,35 +171,54 @@ services::Status PCACorrelationBase<algorithmFPType, cpu>::computeCorrelationEig
 
     DAAL_OVERFLOW_CHECK_BY_MULTIPLICATION(size_t, nFeatures, nFeatures);
     DAAL_OVERFLOW_CHECK_BY_MULTIPLICATION(size_t, nFeatures * nFeatures, sizeof(algorithmFPType));
+    TArray<algorithmFPType, cpu> matrixCopy(nFeatures * nFeatures);
+    DAAL_CHECK_MALLOC(matrixCopy.get());
 
-    TArray<algorithmFPType, cpu> fullEigenvectors(nFeatures * nFeatures);
-    DAAL_CHECK_MALLOC(fullEigenvectors.get());
-    algorithmFPType * fullEigenvectorsArray = fullEigenvectors.get();
-
-    TArray<algorithmFPType, cpu> fullEigenvalues(nFeatures);
-    DAAL_CHECK_MALLOC(fullEigenvalues.get());
-    algorithmFPType * fullEigenvaluesArray = fullEigenvalues.get();
-
-    copyArray(nFeatures * nFeatures, correlationArray, fullEigenvectorsArray);
-
-    services::Status s = computeEigenvectorsInplace(nFeatures, fullEigenvectorsArray, fullEigenvaluesArray);
-    DAAL_CHECK_STATUS_VAR(s);
-
-    s = sortEigenvectorsDescending(nFeatures, fullEigenvectorsArray, fullEigenvaluesArray);
-    DAAL_CHECK_STATUS_VAR(s);
+    algorithmFPType * matrixArray = matrixCopy.get();
+    copyArray(nFeatures * nFeatures, correlationArray, matrixArray);
 
     WriteOnlyRows<algorithmFPType, cpu> eigenvectorsBlock(eigenvectors, 0, nComponents);
     DAAL_CHECK_BLOCK_STATUS(eigenvectorsBlock);
     algorithmFPType * eigenvectorsArray = eigenvectorsBlock.get();
 
+    TArray<algorithmFPType, cpu> fullEigenvalues(nFeatures);
+    DAAL_CHECK_MALLOC(fullEigenvalues.get());
+    algorithmFPType * fullEigenvaluesArray = fullEigenvalues.get();
+
     WriteOnlyRows<algorithmFPType, cpu> eigenvaluesBlock(eigenvalues, 0, 1);
     DAAL_CHECK_BLOCK_STATUS(eigenvaluesBlock);
     algorithmFPType * eigenvaluesArray = eigenvaluesBlock.get();
 
-    copyArray(nFeatures * nComponents, fullEigenvectorsArray, eigenvectorsArray);
-    copyArray(nComponents, fullEigenvaluesArray, eigenvaluesArray);
+    // SYEVR branch
+    // In this case, we compute only nComponents eigenvectors and then sort them in descending order
+    // inside the 'computeEigenvectorsInplaceSyevr' function
+    if (nComponents < 0.25 * nFeatures)
+    {
+        services::Status s = computeEigenvectorsInplaceSyevr(nFeatures, nComponents, matrixArray, fullEigenvaluesArray);
+        DAAL_CHECK_STATUS_VAR(s);
+
+        s = sortEigenvectorsDescending(nFeatures, nComponents, matrixArray, fullEigenvaluesArray);
+        DAAL_CHECK_STATUS_VAR(s);
+
+        copyArray(nFeatures * nComponents, matrixArray, eigenvectorsArray);
+        copyArray(nComponents, fullEigenvaluesArray, eigenvaluesArray);
+        return s;
+    }
+    // SYEVD branch
+    // In this case, we compute all eigenvectors and then sort all of them in descending order
+    // and copy the first nComponents eigenvectors to the output
+    else
+    {
+        services::Status s = computeEigenvectorsInplace(nFeatures, matrixArray, fullEigenvaluesArray);
+        DAAL_CHECK_STATUS_VAR(s);
+
+        s = sortEigenvectorsDescending(nFeatures, nFeatures, matrixArray, fullEigenvaluesArray);
+        DAAL_CHECK_STATUS_VAR(s);
 
-    return s;
+        copyArray(nFeatures * nComponents, matrixArray, eigenvectorsArray);
+        copyArray(nComponents, fullEigenvaluesArray, eigenvaluesArray);
+        return s;
+    }
 }
 
 template <typename algorithmFPType, CpuType cpu>
@@ -223,6 +242,59 @@ services::Status PCACorrelationBase<algorithmFPType, cpu>::computeEigenvectorsIn
     return services::Status();
 }
 
+template <typename algorithmFPType, CpuType cpu>
+services::Status PCACorrelationBase<algorithmFPType, cpu>::computeEigenvectorsInplaceSyevr(size_t nFeatures, size_t nComponents,
+                                                                                           algorithmFPType * eigenvectors,
+                                                                                           algorithmFPType * eigenvalues)
+{
+    char jobz  = 'V';
+    char range = 'I';
+    char uplo  = 'U';
+
+    DAAL_INT il = static_cast<DAAL_INT>(nFeatures - nComponents + 1);
+    DAAL_INT iu = static_cast<DAAL_INT>(nFeatures);
+    DAAL_INT m;
+    DAAL_INT info;
+
+    algorithmFPType abstol = -1;
+
+    // Workspace query
+    algorithmFPType work_query;
+    DAAL_INT iwork_query;
+    DAAL_INT lwork  = -1;
+    DAAL_INT liwork = -1;
+
+    TArray<DAAL_INT, cpu> dummy_isuppz(2);
+    DAAL_CHECK_MALLOC(dummy_isuppz.get());
+
+    LapackInst<algorithmFPType, cpu>::xsyevr(&jobz, &range, &uplo, (DAAL_INT *)(&nFeatures), eigenvectors, (DAAL_INT *)(&nFeatures), nullptr, nullptr,
+                                             &il, &iu, &abstol, &m, nullptr, nullptr, (DAAL_INT *)(&nFeatures), dummy_isuppz.get(), &work_query,
+                                             &lwork, &iwork_query, &liwork, &info);
+
+    lwork  = static_cast<DAAL_INT>(work_query);
+    liwork = iwork_query;
+
+    TArray<algorithmFPType, cpu> work(lwork);
+    TArray<DAAL_INT, cpu> iwork(liwork);
+    TArray<DAAL_INT, cpu> isuppz(2 * nComponents);
+    DAAL_CHECK_MALLOC(work.get() && iwork.get() && isuppz.get());
+
+    TArray<algorithmFPType, cpu> temp_eigenvectors(nFeatures * nFeatures);
+
+    LapackInst<algorithmFPType, cpu>::xsyevr(&jobz, &range, &uplo, (DAAL_INT *)(&nFeatures), eigenvectors, (DAAL_INT *)(&nFeatures), nullptr, nullptr,
+                                             &il, &iu, &abstol, &m, eigenvalues, temp_eigenvectors.get(), (DAAL_INT *)(&nFeatures), isuppz.get(),
+                                             work.get(), &lwork, iwork.get(), &liwork, &info);
+
+    if (info != 0 || m != static_cast<DAAL_INT>(nComponents))
+    {
+        return services::Status(services::ErrorPCAFailedToComputeCorrelationEigenvalues);
+    }
+
+    copyArray(nFeatures * nComponents, temp_eigenvectors.get(), eigenvectors);
+
+    return services::Status();
+}
+
 template <typename algorithmFPType, CpuType cpu>
 services::Status PCACorrelationBase<algorithmFPType, cpu>::sortEigenvectorsDescending(size_t nFeatures, algorithmFPType * eigenvectors,
                                                                                       algorithmFPType * eigenvalues)
@@ -245,6 +317,46 @@ services::Status PCACorrelationBase<algorithmFPType, cpu>::sortEigenvectorsDesce
     return services::Status();
 }
 
+template <typename algorithmFPType, CpuType cpu>
+services::Status PCACorrelationBase<algorithmFPType, cpu>::sortEigenvectorsDescending(size_t nFeatures, size_t nComponents,
+                                                                                      algorithmFPType * eigenvectors, algorithmFPType * eigenvalues)
+{
+    const algorithmFPType eps = std::numeric_limits<algorithmFPType>::epsilon();
+
+    for (size_t i = 0; i < nComponents / 2; i++)
+    {
+        algorithmFPType tmp     = eigenvalues[i];
+        algorithmFPType tmp_rev = eigenvalues[nComponents - 1 - i];
+
+        if (tmp < eps) tmp = algorithmFPType(0);
+        if (tmp_rev < eps) tmp_rev = algorithmFPType(0);
+
+        eigenvalues[i]                   = tmp_rev;
+        eigenvalues[nComponents - 1 - i] = tmp;
+    }
+
+    if (nComponents % 2 != 0)
+    {
+        size_t mid = nComponents / 2;
+        if (eigenvalues[mid] < eps)
+        {
+            eigenvalues[mid] = algorithmFPType(0);
+        }
+    }
+
+    TArray<algorithmFPType, cpu> eigenvectorTmp(nFeatures);
+    DAAL_CHECK_MALLOC(eigenvectorTmp.get());
+
+    for (size_t i = 0; i < nComponents / 2; i++)
+    {
+        copyArray(nFeatures, eigenvectors + i * nFeatures, eigenvectorTmp.get());
+        copyArray(nFeatures, eigenvectors + nFeatures * (nComponents - 1 - i), eigenvectors + i * nFeatures);
+        copyArray(nFeatures, eigenvectorTmp.get(), eigenvectors + nFeatures * (nComponents - 1 - i));
+    }
+
+    return services::Status();
+}
+
 } // namespace internal
 } // namespace pca
 } // namespace algorithms

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

@@ -163,6 +163,17 @@ static train_result<Task> call_daal_kernel_finalize_train(const context_cpu& ctx
         result.set_variances(homogen_table::wrap(arr_vars, 1, column_count));
     }
 
+    if (desc.get_result_options().test(result_options::noise_variance)) {
+        double noiseVariance = 0.0;
+        interop::status_to_exception(dal::backend::dispatch_by_cpu(ctx, [&](auto cpu) {
+            return daal_pca_cor_kernel_t<
+                       Float,
+                       dal::backend::interop::to_daal_cpu_type<decltype(cpu)>::value>()
+                .computeNoiseVariances(*daal_eigenvalues, *daal_variances, noiseVariance);
+        }));
+        result.set_noise_variance(noiseVariance);
+    }
+
     if (desc.get_result_options().test(result_options::explained_variances_ratio)) {
         result.set_explained_variances_ratio(
             homogen_table::wrap(arr_explained_variances_ratio, 1, component_count));

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

@@ -126,6 +126,17 @@ static result_t call_daal_kernel(const context_cpu& ctx,
         result.set_singular_values(homogen_table::wrap(arr_singular_values, 1, component_count));
     }
 
+    if (desc.get_result_options().test(result_options::noise_variance)) {
+        double noiseVariance = 0.0;
+        interop::status_to_exception(dal::backend::dispatch_by_cpu(ctx, [&](auto cpu) {
+            return daal_pca_cor_kernel_t<
+                       Float,
+                       dal::backend::interop::to_daal_cpu_type<decltype(cpu)>::value>()
+                .computeNoiseVariances(*daal_eigenvalues, *daal_variances, noiseVariance);
+        }));
+        result.set_noise_variance(noiseVariance);
+    }
+
     if (desc.get_result_options().test(result_options::explained_variances_ratio)) {
         result.set_explained_variances_ratio(
             homogen_table::wrap(arr_explained_variances_ratio, 1, component_count));

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

@@ -98,6 +98,17 @@ static result_t call_daal_kernel(const context_cpu& ctx,
         result.set_eigenvectors(homogen_table::wrap(arr_eigvec, component_count, column_count));
     }
 
+    if (desc.get_result_options().test(result_options::noise_variance)) {
+        double noiseVariance = 0.0;
+        interop::status_to_exception(dal::backend::dispatch_by_cpu(ctx, [&](auto cpu) {
+            return daal_pca_cor_kernel_t<
+                       Float,
+                       dal::backend::interop::to_daal_cpu_type<decltype(cpu)>::value>()
+                .computeNoiseVariances(*daal_eigenvalues, *daal_variances, noiseVariance);
+        }));
+        result.set_noise_variance(noiseVariance);
+    }
+
     if (desc.get_result_options().test(result_options::eigenvalues)) {
         result.set_eigenvalues(homogen_table::wrap(arr_eigval, 1, component_count));
     }

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

@@ -142,6 +142,17 @@ static result_t call_daal_kernel(const context_cpu& ctx,
         result.set_variances(homogen_table::wrap(arr_vars, 1, column_count));
     }
 
+    if (desc.get_result_options().test(result_options::noise_variance)) {
+        double noiseVariance = 0.0;
+        interop::status_to_exception(dal::backend::dispatch_by_cpu(ctx, [&](auto cpu) {
+            return daal_pca_svd_kernel_t<
+                       Float,
+                       dal::backend::interop::to_daal_cpu_type<decltype(cpu)>::value>()
+                .computeNoiseVariances(*daal_eigenvalues, *daal_variances, noiseVariance);
+        }));
+        result.set_noise_variance(noiseVariance);
+    }
+
     if (desc.get_result_options().test(result_options::means)) {
         result.set_means(homogen_table::wrap(arr_means, 1, column_count));
     }