ENH: Non-batched linear regression for high-dimensional problems

cpp/daal/src/algorithms/linear_model/linear_model_hyperparameter_impl.h

@@ -41,8 +41,11 @@ namespace internal
  */
 enum HyperparameterId
 {
-    denseUpdateStepBlockSize = 0,
-    hyperparameterIdCount    = denseUpdateStepBlockSize + 1
+    denseUpdateStepBlockSize     = 0,
+    denseUpdateMaxColsBatched    = 1,
+    denseSmallRowsThreshold      = 2,
+    denseSmallRowsMaxColsBatched = 3,
+    hyperparameterIdCount        = 4
 };
 
 enum DoubleHyperparameterId

cpp/daal/src/algorithms/linear_model/linear_model_train_normeq_update_impl.i

@@ -173,6 +173,59 @@ Status UpdateKernel<algorithmFPType, cpu>::compute(const NumericTable & xTable,
     return UpdateKernel<algorithmFPType, cpu>::compute(xTable, yTable, xtxTable, xtyTable, initializeResult, interceptFlag, nullptr);
 }
 
+template <typename algorithmFPType, CpuType cpu>
+Status computeNonBatchedAggregates(const DAAL_INT nRows, const DAAL_INT nCols, const DAAL_INT nResponses, bool initializeResult, bool interceptFlag,
+                                   const algorithmFPType * xPtr, const algorithmFPType * yPtr, algorithmFPType * xtx, algorithmFPType * xty)
+{
+    DAAL_INT nBetasIntercept            = nCols + static_cast<int>(interceptFlag);
+    DAAL_INT one_int                    = 1;
+    algorithmFPType one                 = 1;
+    algorithmFPType zero                = 0;
+    algorithmFPType * mult_current_data = initializeResult ? &zero : &one;
+    daal::services::internal::TArray<algorithmFPType, cpu> ones;
+    if (interceptFlag)
+    {
+        ones = daal::services::internal::TArray<algorithmFPType, cpu>(nRows);
+        std::fill(ones.get(), ones.get() + nRows, algorithmFPType(1));
+    }
+
+    BlasInst<algorithmFPType, cpu>::xsyrk("U", "N", &nCols, &nRows, &one, xPtr, &nCols, mult_current_data, xtx, &nBetasIntercept);
+    if (interceptFlag)
+    {
+        BlasInst<algorithmFPType, cpu>::xgemv("N", &nCols, &nRows, &one, xPtr, &nCols, ones.get(), &one_int, mult_current_data,
+                                              xtx + static_cast<size_t>(nBetasIntercept) * static_cast<size_t>(nCols), &one_int);
+        const size_t idx_last = static_cast<size_t>(nBetasIntercept) * static_cast<size_t>(nBetasIntercept) - 1;
+        if (initializeResult)
+            xtx[idx_last] = nRows;
+        else
+            xtx[idx_last] += nRows;
+    }
+
+    if (nResponses == 1)
+    {
+        BlasInst<algorithmFPType, cpu>::xgemv("N", &nCols, &nRows, &one, xPtr, &nCols, yPtr, &one_int, mult_current_data, xty, &one_int);
+        if (interceptFlag)
+        {
+            const algorithmFPType last_val = BlasInst<algorithmFPType, cpu>::xxdot(&nRows, yPtr, &one_int, ones.get(), &one_int);
+            if (initializeResult)
+                xty[nCols] = last_val;
+            else
+                xty[nCols] += last_val;
+        }
+    }
+    else
+    {
+        BlasInst<algorithmFPType, cpu>::xgemm("N", "T", &nCols, &nResponses, &nRows, &one, xPtr, &nCols, yPtr, &nResponses, mult_current_data, xty,
+                                              &nBetasIntercept);
+        if (interceptFlag)
+        {
+            BlasInst<algorithmFPType, cpu>::xgemv("N", &nResponses, &nRows, &one, yPtr, &nResponses, ones.get(), &one_int, mult_current_data,
+                                                  xty + nCols, &nBetasIntercept);
+        }
+    }
+    return Status();
+}
+
 template <typename algorithmFPType, CpuType cpu>
 Status UpdateKernel<algorithmFPType, cpu>::compute(const NumericTable & xTable, const NumericTable & yTable, NumericTable & xtxTable,
                                                    NumericTable & xtyTable, bool initializeResult, bool interceptFlag,
@@ -193,6 +246,47 @@ Status UpdateKernel<algorithmFPType, cpu>::compute(const NumericTable & xTable,
     DAAL_CHECK_BLOCK_STATUS(xtyBlock);
     algorithmFPType * xty = xtyBlock.get();
 
+    /// Logic here is as follows: it needs to compute t(X)*X and t(X)*y.
+    /// If both are done together, it's possible to reuse caches of data to speed up computations,
+    /// which the code here does by dividing the data into batches of rows on which both aggregates
+    /// are computed, with the batches processed in parallel. But as the number of columns in the
+    /// data grows, the potential speed gains from calculating both aggregates simultaneously
+    /// decreases, and the memory requirements increase, which can become a problem when there are
+    /// many threads in the system. Hence, if the number of columns is too large, it will compute
+    /// both aggregates independently, in separate calls to BLAS functions, while if the number of
+    /// columns is reasonably small, will prefer the batched procedure which typically ends up
+    /// being faster.
+    DAAL_INT64 maxColsBatched          = 4096;
+    DAAL_INT64 smallRowsThreshold      = 10000;
+    DAAL_INT64 smallRowsMaxColsBatched = 1024;
+    if (hyperparameter != nullptr)
+    {
+        services::Status status = hyperparameter->find(denseUpdateMaxColsBatched, maxColsBatched);
+        DAAL_CHECK(maxColsBatched > 0, services::ErrorIncorrectDataRange);
+        DAAL_CHECK_STATUS_VAR(status);
+
+        status = hyperparameter->find(denseSmallRowsThreshold, smallRowsThreshold);
+        DAAL_CHECK(smallRowsThreshold > 0, services::ErrorIncorrectDataRange);
+        DAAL_CHECK_STATUS_VAR(status);
+
+        status = hyperparameter->find(denseSmallRowsMaxColsBatched, smallRowsMaxColsBatched);
+        DAAL_CHECK(smallRowsMaxColsBatched > 0, services::ErrorIncorrectDataRange);
+        DAAL_CHECK_STATUS_VAR(status);
+    }
+
+    const bool use_non_batched_route = nBetas >= maxColsBatched || (nRows >= smallRowsThreshold && nBetas >= smallRowsMaxColsBatched);
+    if (use_non_batched_route)
+    {
+        const DAAL_INT nCols = xTable.getNumberOfColumns();
+        ReadRowsType xBlock(const_cast<NumericTable &>(xTable), 0, nRows);
+        DAAL_CHECK_BLOCK_STATUS(xBlock);
+        ReadRowsType yBlock(const_cast<NumericTable &>(yTable), 0, nRows);
+        DAAL_CHECK_BLOCK_STATUS(yBlock);
+        const algorithmFPType * xPtr = xBlock.get();
+        const algorithmFPType * yPtr = yBlock.get();
+        return computeNonBatchedAggregates<algorithmFPType, cpu>(nRows, nCols, nResponses, initializeResult, interceptFlag, xPtr, yPtr, xtx, xty);
+    }
+
     /* Initialize output arrays by zero in case of batch mode */
     if (initializeResult)
     {

cpp/daal/src/algorithms/linear_regression/linear_regression_hyperparameter.cpp

@@ -64,17 +64,40 @@ services::Status convert(const Hyperparameter * params, services::SharedPtr<Line
 
     if (params != nullptr)
     {
-        DAAL_INT64 denseUpdateStepBlockSize = 0l;
+        DAAL_INT64 denseUpdateStepBlockSize     = 0l;
+        DAAL_INT64 denseUpdateMaxColsBatched    = 0l;
+        DAAL_INT64 denseSmallRowsThreshold      = 0l;
+        DAAL_INT64 denseSmallRowsMaxColsBatched = 0l;
 
         auto * const resultPtr = new LinearModelHyperparameter();
         DAAL_CHECK_MALLOC(resultPtr);
         result.reset(resultPtr);
 
+        /// Getters
         st |= params->find(HyperparameterId::denseUpdateStepBlockSize, denseUpdateStepBlockSize);
         DAAL_CHECK(st, services::ErrorHyperparameterNotFound);
 
+        st |= params->find(HyperparameterId::denseUpdateMaxColsBatched, denseUpdateMaxColsBatched);
+        DAAL_CHECK(st, services::ErrorHyperparameterNotFound);
+
+        st |= params->find(HyperparameterId::denseSmallRowsThreshold, denseSmallRowsThreshold);
+        DAAL_CHECK(st, services::ErrorHyperparameterNotFound);
+
+        st |= params->find(HyperparameterId::denseSmallRowsMaxColsBatched, denseSmallRowsMaxColsBatched);
+        DAAL_CHECK(st, services::ErrorHyperparameterNotFound);
+
+        /// Setters
         st |= result->set(LinearModelHyperparameterId::denseUpdateStepBlockSize, denseUpdateStepBlockSize);
         DAAL_CHECK(st, services::ErrorHyperparameterCanNotBeSet);
+
+        st |= result->set(LinearModelHyperparameterId::denseUpdateMaxColsBatched, denseUpdateMaxColsBatched);
+        DAAL_CHECK(st, services::ErrorHyperparameterCanNotBeSet);
+
+        st |= result->set(LinearModelHyperparameterId::denseSmallRowsThreshold, denseSmallRowsThreshold);
+        DAAL_CHECK(st, services::ErrorHyperparameterCanNotBeSet);
+
+        st |= result->set(LinearModelHyperparameterId::denseSmallRowsMaxColsBatched, denseSmallRowsMaxColsBatched);
+        DAAL_CHECK(st, services::ErrorHyperparameterCanNotBeSet);
     }
     else
     {

cpp/daal/src/algorithms/linear_regression/linear_regression_hyperparameter_impl.h

@@ -45,8 +45,11 @@ typedef linear_model::internal::Hyperparameter LinearModelHyperparameter;
  */
 enum HyperparameterId
 {
-    denseUpdateStepBlockSize = 0,
-    hyperparameterIdCount    = denseUpdateStepBlockSize + 1
+    denseUpdateStepBlockSize     = 0,
+    denseUpdateMaxColsBatched    = 1,
+    denseSmallRowsThreshold      = 2,
+    denseSmallRowsMaxColsBatched = 3,
+    hyperparameterIdCount        = 4
 };
 
 enum DoubleHyperparameterId

cpp/daal/src/externals/service_blas.h

@@ -103,14 +103,14 @@ struct Blas
 {
     typedef typename _impl<fpType, cpu>::SizeType SizeType;
 
-    static void xsyrk(char * uplo, char * trans, SizeType * p, SizeType * n, fpType * alpha, fpType * a, SizeType * lda, fpType * beta, fpType * ata,
-                      SizeType * ldata)
+    static void xsyrk(const char * uplo, const char * trans, const SizeType * p, const SizeType * n, const fpType * alpha, const fpType * a,
+                      const SizeType * lda, const fpType * beta, fpType * ata, const SizeType * ldata)
     {
         _impl<fpType, cpu>::xsyrk(uplo, trans, p, n, alpha, a, lda, beta, ata, ldata);
     }
 
-    static void xxsyrk(char * uplo, char * trans, SizeType * p, SizeType * n, fpType * alpha, fpType * a, SizeType * lda, fpType * beta, fpType * ata,
-                       SizeType * ldata)
+    static void xxsyrk(const char * uplo, const char * trans, const SizeType * p, const SizeType * n, const fpType * alpha, const fpType * a,
+                       const SizeType * lda, const fpType * beta, fpType * ata, const SizeType * ldata)
     {
         _impl<fpType, cpu>::xxsyrk(uplo, trans, p, n, alpha, a, lda, beta, ata, ldata);
     }

cpp/daal/src/externals/service_blas_declar_ref.h

@@ -34,10 +34,10 @@ namespace ref
 {
 extern "C"
 {
-    extern void ssyrk_(const char *, const char *, const DAAL_INT *, const DAAL_INT *, const float *, float *, const DAAL_INT *, const float *,
+    extern void ssyrk_(const char *, const char *, const DAAL_INT *, const DAAL_INT *, const float *, const float *, const DAAL_INT *, const float *,
                        float *, const DAAL_INT *);
-    extern void dsyrk_(const char *, const char *, const DAAL_INT *, const DAAL_INT *, const double *, double *, const DAAL_INT *, const double *,
-                       double *, const DAAL_INT *);
+    extern void dsyrk_(const char *, const char *, const DAAL_INT *, const DAAL_INT *, const double *, const double *, const DAAL_INT *,
+                       const double *, double *, const DAAL_INT *);
 
     extern void ssyr_(const char *, const DAAL_INT *, const float *, const float *, const DAAL_INT *, float *, const DAAL_INT *);
     extern void dsyr_(const char *, const DAAL_INT *, const double *, const double *, const DAAL_INT *, double *, const DAAL_INT *);

cpp/daal/src/externals/service_blas_mkl.h

@@ -50,14 +50,14 @@ struct MklBlas<double, cpu>
 {
     typedef DAAL_INT SizeType;
 
-    static void xsyrk(char * uplo, char * trans, DAAL_INT * p, DAAL_INT * n, double * alpha, double * a, DAAL_INT * lda, double * beta, double * ata,
-                      DAAL_INT * ldata)
+    static void xsyrk(const char * uplo, const char * trans, const DAAL_INT * p, const DAAL_INT * n, const double * alpha, const double * a,
+                      const DAAL_INT * lda, const double * beta, double * ata, const DAAL_INT * ldata)
     {
         __DAAL_MKLFN_CALL_BLAS(dsyrk, (uplo, trans, (MKL_INT *)p, (MKL_INT *)n, alpha, a, (MKL_INT *)lda, beta, ata, (MKL_INT *)ldata));
     }
 
-    static void xxsyrk(char * uplo, char * trans, DAAL_INT * p, DAAL_INT * n, double * alpha, double * a, DAAL_INT * lda, double * beta, double * ata,
-                       DAAL_INT * ldata)
+    static void xxsyrk(const char * uplo, const char * trans, const DAAL_INT * p, const DAAL_INT * n, const double * alpha, const double * a,
+                       const DAAL_INT * lda, const double * beta, double * ata, const DAAL_INT * ldata)
     {
         int old_nthr = mkl_set_num_threads_local(1);
         __DAAL_MKLFN_CALL_BLAS(dsyrk, (uplo, trans, (MKL_INT *)p, (MKL_INT *)n, alpha, a, (MKL_INT *)lda, beta, ata, (MKL_INT *)ldata));
@@ -155,14 +155,14 @@ struct MklBlas<float, cpu>
 {
     typedef DAAL_INT SizeType;
 
-    static void xsyrk(char * uplo, char * trans, DAAL_INT * p, DAAL_INT * n, float * alpha, float * a, DAAL_INT * lda, float * beta, float * ata,
-                      DAAL_INT * ldata)
+    static void xsyrk(const char * uplo, const char * trans, const DAAL_INT * p, const DAAL_INT * n, const float * alpha, const float * a,
+                      const DAAL_INT * lda, const float * beta, float * ata, const DAAL_INT * ldata)
     {
         __DAAL_MKLFN_CALL_BLAS(ssyrk, (uplo, trans, (MKL_INT *)p, (MKL_INT *)n, alpha, a, (MKL_INT *)lda, beta, ata, (MKL_INT *)ldata));
     }
 
-    static void xxsyrk(char * uplo, char * trans, DAAL_INT * p, DAAL_INT * n, float * alpha, float * a, DAAL_INT * lda, float * beta, float * ata,
-                       DAAL_INT * ldata)
+    static void xxsyrk(const char * uplo, const char * trans, const DAAL_INT * p, const DAAL_INT * n, const float * alpha, const float * a,
+                       const DAAL_INT * lda, const float * beta, float * ata, const DAAL_INT * ldata)
     {
         int old_nthr = mkl_set_num_threads_local(1);
         __DAAL_MKLFN_CALL_BLAS(ssyrk, (uplo, trans, (MKL_INT *)p, (MKL_INT *)n, alpha, a, (MKL_INT *)lda, beta, ata, (MKL_INT *)ldata));

cpp/daal/src/externals/service_blas_ref.h

@@ -46,14 +46,14 @@ struct OpenBlas<double, cpu>
 {
     typedef DAAL_INT SizeType;
 
-    static void xsyrk(char * uplo, char * trans, DAAL_INT * p, DAAL_INT * n, double * alpha, double * a, DAAL_INT * lda, double * beta, double * ata,
-                      DAAL_INT * ldata)
+    static void xsyrk(const char * uplo, const char * trans, const DAAL_INT * p, const DAAL_INT * n, const double * alpha, const double * a,
+                      const DAAL_INT * lda, const double * beta, double * ata, const DAAL_INT * ldata)
     {
         dsyrk_(uplo, trans, p, n, alpha, a, lda, beta, ata, ldata);
     }
 
-    static void xxsyrk(char * uplo, char * trans, DAAL_INT * p, DAAL_INT * n, double * alpha, double * a, DAAL_INT * lda, double * beta, double * ata,
-                       DAAL_INT * ldata)
+    static void xxsyrk(const char * uplo, const char * trans, const DAAL_INT * p, const DAAL_INT * n, const double * alpha, const double * a,
+                       const DAAL_INT * lda, const double * beta, double * ata, const DAAL_INT * ldata)
     {
         openblas_thread_setter ots(1);
         dsyrk_(uplo, trans, p, n, alpha, a, lda, beta, ata, ldata);
@@ -138,14 +138,14 @@ struct OpenBlas<float, cpu>
 {
     typedef DAAL_INT SizeType;
 
-    static void xsyrk(char * uplo, char * trans, DAAL_INT * p, DAAL_INT * n, float * alpha, float * a, DAAL_INT * lda, float * beta, float * ata,
-                      DAAL_INT * ldata)
+    static void xsyrk(const char * uplo, const char * trans, const DAAL_INT * p, const DAAL_INT * n, const float * alpha, const float * a,
+                      const DAAL_INT * lda, const float * beta, float * ata, const DAAL_INT * ldata)
     {
         ssyrk_(uplo, trans, p, n, alpha, a, lda, beta, ata, ldata);
     }
 
-    static void xxsyrk(char * uplo, char * trans, DAAL_INT * p, DAAL_INT * n, float * alpha, float * a, DAAL_INT * lda, float * beta, float * ata,
-                       DAAL_INT * ldata)
+    static void xxsyrk(const char * uplo, const char * trans, const DAAL_INT * p, const DAAL_INT * n, const float * alpha, const float * a,
+                       const DAAL_INT * lda, const float * beta, float * ata, const DAAL_INT * ldata)
     {
         openblas_thread_setter ots(1);
         ssyrk_(uplo, trans, p, n, alpha, a, lda, beta, ata, ldata);

cpp/oneapi/dal/algo/linear_regression/backend/cpu/finalize_train_kernel_norm_eq.cpp

@@ -55,10 +55,21 @@ static daal_lr_hyperparameters_t convert_parameters(const detail::train_paramete
     using daal_lr::internal::HyperparameterId;
 
     const std::int64_t block = params.get_cpu_macro_block();
+    const std::int64_t max_cols_batched = params.get_cpu_max_cols_batched();
+    const std::int64_t small_rows_threshold = params.get_cpu_small_rows_threshold();
+    const std::int64_t small_rows_max_cols_batched = params.get_cpu_small_rows_max_cols_batched();
 
     daal_lr_hyperparameters_t daal_hyperparameter;
     auto status = daal_hyperparameter.set(HyperparameterId::denseUpdateStepBlockSize, block);
     interop::status_to_exception(status);
+    status = daal_hyperparameter.set(HyperparameterId::denseUpdateMaxColsBatched, max_cols_batched);
+    interop::status_to_exception(status);
+    status =
+        daal_hyperparameter.set(HyperparameterId::denseSmallRowsThreshold, small_rows_threshold);
+    interop::status_to_exception(status);
+    status = daal_hyperparameter.set(HyperparameterId::denseSmallRowsMaxColsBatched,
+                                     small_rows_max_cols_batched);
+    interop::status_to_exception(status);
 
     return daal_hyperparameter;
 }

cpp/oneapi/dal/algo/linear_regression/backend/cpu/partial_train_kernel_norm_eq.cpp

@@ -47,10 +47,21 @@ static daal_hyperparameters_t convert_parameters(const detail::train_parameters<
     using daal_lr::internal::HyperparameterId;
 
     const std::int64_t block = params.get_cpu_macro_block();
+    const std::int64_t max_cols_batched = params.get_cpu_max_cols_batched();
+    const std::int64_t small_rows_threshold = params.get_cpu_small_rows_threshold();
+    const std::int64_t small_rows_max_cols_batched = params.get_cpu_small_rows_max_cols_batched();
 
     daal_hyperparameters_t daal_hyperparameter;
     auto status = daal_hyperparameter.set(HyperparameterId::denseUpdateStepBlockSize, block);
     interop::status_to_exception(status);
+    status = daal_hyperparameter.set(HyperparameterId::denseUpdateMaxColsBatched, max_cols_batched);
+    interop::status_to_exception(status);
+    status =
+        daal_hyperparameter.set(HyperparameterId::denseSmallRowsThreshold, small_rows_threshold);
+    interop::status_to_exception(status);
+    status = daal_hyperparameter.set(HyperparameterId::denseSmallRowsMaxColsBatched,
+                                     small_rows_max_cols_batched);
+    interop::status_to_exception(status);
 
     return daal_hyperparameter;
 }

cpp/oneapi/dal/algo/linear_regression/backend/cpu/train_kernel_norm_eq.cpp

@@ -58,10 +58,21 @@ static daal_lr_hyperparameters_t convert_parameters(const detail::train_paramete
     using daal_lr::internal::HyperparameterId;
 
     const std::int64_t block = params.get_cpu_macro_block();
+    const std::int64_t max_cols_batched = params.get_cpu_max_cols_batched();
+    const std::int64_t small_rows_threshold = params.get_cpu_small_rows_threshold();
+    const std::int64_t small_rows_max_cols_batched = params.get_cpu_small_rows_max_cols_batched();
 
     daal_lr_hyperparameters_t daal_hyperparameter;
     auto status = daal_hyperparameter.set(HyperparameterId::denseUpdateStepBlockSize, block);
     interop::status_to_exception(status);
+    status = daal_hyperparameter.set(HyperparameterId::denseUpdateMaxColsBatched, max_cols_batched);
+    interop::status_to_exception(status);
+    status =
+        daal_hyperparameter.set(HyperparameterId::denseSmallRowsThreshold, small_rows_threshold);
+    interop::status_to_exception(status);
+    status = daal_hyperparameter.set(HyperparameterId::denseSmallRowsMaxColsBatched,
+                                     small_rows_max_cols_batched);
+    interop::status_to_exception(status);
 
     return daal_hyperparameter;
 }

cpp/oneapi/dal/algo/linear_regression/test/batch.cpp

@@ -39,6 +39,9 @@ class lr_batch_test : public lr_test<TestType, lr_batch_test<TestType>> {
         this->f_count_ = GENERATE(2, 17);
         this->r_count_ = GENERATE(2, 15);
         this->intercept_ = GENERATE(0, 1);
+        if (this->get_policy().is_cpu()) {
+            this->use_non_batched_route = GENERATE(false, true);
+        }
     }
 };