Полностью развернутый QR быстрее, чем неповоротный QR в С ++? (с библиотекой Eigen) - PullRequest
Новая
       85

Полностью развернутый QR быстрее, чем неповоротный QR в С ++? (с библиотекой Eigen)

0 голосов
/ 26 мая 2020

Я получаю странные результаты производительности при сравнении различных реализаций QR Eigen на C ++ [Я обращаюсь к среде C ++ с помощью пакета RcppEigen из R - но строка sr c - это просто код C ++]: 

# 1) Non-pivoted QR decomposition:

src<-'
using Eigen::Map;
using Eigen::MatrixXd;
using Eigen::HouseholderQR;
const Map<MatrixXd> X(as<Map<MatrixXd> >(AA)); // mapping AA to X
int m(X.rows());                                 
MatrixXd Q(m,m);
MatrixXd R(X);)
HouseholderQR<MatrixXd> qr(X);  
Q = qr.householderQ();
R = qr.matrixQR().triangularView<Eigen::Upper>();    
return List::create(Named("Q") = Q,
                    Named("R") = R);
'

HHqr<-cxxfunction(signature(AA="matrix"),body=src,plugin = "RcppEigen")

# 2) Column-pivoted QR decomposition:

src<-'
using Eigen::Map;
using Eigen::MatrixXd;
using Eigen::ColPivHouseholderQR;
const Map<MatrixXd> X(as<Map<MatrixXd> >(AA)); // mapping AA to X 
int m(X.rows());   
int n(X.cols());           
MatrixXd Q(m,m),P(n,n); 
MatrixXd R(X);
ColPivHouseholderQR<MatrixXd> qr(X);                
P = qr.colsPermutation();
Q = qr.householderQ();
R = qr.matrixQR().triangularView<Eigen::Upper>();
return List::create(Named("Q") = Q,
                    Named("R") = R,
                    Named("P") = P,
                    Named("Rank") = qr.rank());
'
CPHHqr<-cxxfunction(signature(AA="matrix"),body=src,plugin = "RcppEigen")

# 3) Full-pivoted QR decomposition:

src<-'
using Eigen::Map;
using Eigen::MatrixXd;
using Eigen::FullPivHouseholderQR;
const Map<MatrixXd> X(as<Map<MatrixXd> >(AA)); // mapping AA to X
m(X.rows());                                 
int n(X.cols());                                         
MatrixXd Q(m,m),P(n,n);                        
MatrixXd R(X);                                 
FullPivHouseholderQR<MatrixXd> qr(X);          
P = qr.colsPermutation();
Q = qr.matrixQ();
R = qr.matrixQR().triangularView<Eigen::Upper>();
return List::create(Named("Q") = Q,
                    Named("R") = R,
                    Named("P") = P,
                    Named("Rank") = qr.rank());
'
FPHHqr<-cxxfunction(signature(AA="matrix"),body=src,plugin = "RcppEigen")

и вот результаты теста: 

library(fBasics)
X<-hilbert(500)   # test matrix
gc()
library(rbenchmark)
benchmark(HHqr(X),
          CPHHqr(X),
          FPHHqr(X),
          order = "elapsed")  

        test replications elapsed relative user.self sys.self user.child sys.child
>3 FPHHqr(X)          100    2.10    1.000      1.65     0.42         NA        NA
>1   HHqr(X)          100    5.61    2.671      5.29     0.22         NA        NA
>2 CPHHqr(X)          100    7.47    3.557      6.93     0.47         NA        NA

Как вообще ПОЛНОСТЬЮ вращающийся QR 2.5x быстрее, чем НЕВЕРНЫЙ? Спасибо.

1 Ответ

0 голосов
/ 27 мая 2020

Может ли это быть функцией библиотек BLAS / LAPACK (даже несмотря на то, что Eigen3 делает (делал?) Там некоторые внутренние вычисления)? Я получаю разные результаты (из слегка измененной версии вашего кода. Я вижу 

R> library(fBasics)

R> X <- hilbert(500)   # test matrix

R> library(microbenchmark)

R> microbenchmark(HHqr(X),
+                CPHHqr(X),
+                FPHHqr(X))
Unit: milliseconds
      expr      min       lq    mean   median      uq      max neval cld
   HHqr(X) 31.16688 32.08694 33.6356 32.85707 34.8813  41.4097   100  b 
 CPHHqr(X) 36.99445 37.94523 40.6429 38.79969 41.2618 133.2360   100   c
 FPHHqr(X)  8.71062  9.27885 11.3358  9.58908 11.4974 104.1892   100 a  
R>

Мой код ниже --- и вы можете запустить его вместе с тестом производительности просто через `R cpp :: source Cpp ("nameOfThatFile. cpp").

И ниже мой sessionInfo(). Я использую OpenBLAS, который является многопоточным, и у меня есть шестиядерный компьютер.

Код 

#include <RcppEigen.h>

using namespace Rcpp;

// [[Rcpp::depends(RcppEigen)]]


// 1) Non-pivoted QR decomposition:
// [[Rcpp::export]]
List HHqr(NumericMatrix AA) {
  using Eigen::Map;
  using Eigen::MatrixXd;
  using Eigen::HouseholderQR;
  const Map<MatrixXd> X(as<Map<MatrixXd> >(AA)); // mapping AA to X
  int m(X.rows());
  MatrixXd Q(m,m);
  MatrixXd R(X);
  HouseholderQR<MatrixXd> qr(X);
  Q = qr.householderQ();
  R = qr.matrixQR().triangularView<Eigen::Upper>();
  return List::create(Named("Q") = Q,
                      Named("R") = R);
}

//HHqr<-cxxfunction(signature(AA="matrix"),body=src,plugin = "RcppEigen")

// # 2) Column-pivoted QR decomposition:
// [[Rcpp::export]]
List CPHHqr(NumericMatrix AA) {
  using Eigen::Map;
  using Eigen::MatrixXd;
  using Eigen::ColPivHouseholderQR;
  const Map<MatrixXd> X(as<Map<MatrixXd> >(AA)); // mapping AA to X
  int m(X.rows());
  int n(X.cols());
  MatrixXd Q(m,m),P(n,n);
  MatrixXd R(X);
  ColPivHouseholderQR<MatrixXd> qr(X);
  P = qr.colsPermutation();
  Q = qr.householderQ();
  R = qr.matrixQR().triangularView<Eigen::Upper>();
  return List::create(Named("Q") = Q,
                      Named("R") = R,
                      Named("P") = P,
                      Named("Rank") = qr.rank());
}

//CPHHqr<-cxxfunction(signature(AA="matrix"),body=src,plugin = "RcppEigen")

//# 3) Full-pivoted QR decomposition:
// [[Rcpp::export]]
List FPHHqr(NumericMatrix AA) {
  using Eigen::Map;
  using Eigen::MatrixXd;
  using Eigen::FullPivHouseholderQR;
  const Map<MatrixXd> X(as<Map<MatrixXd> >(AA)); // mapping AA to X
  int m(X.rows());
  int n(X.cols());
  MatrixXd Q(m,m),P(n,n);
  MatrixXd R(X);
  FullPivHouseholderQR<MatrixXd> qr(X);
  P = qr.colsPermutation();
  Q = qr.matrixQ();
  R = qr.matrixQR().triangularView<Eigen::Upper>();
  return List::create(Named("Q") = Q,
                      Named("R") = R,
                      Named("P") = P,
                      Named("Rank") = qr.rank());
}
//FPHHqr<-cxxfunction(signature(AA="matrix"),body=src,plugin = "RcppEigen")

/*** R
library(fBasics)
X <- hilbert(500)   # test matrix
library(microbenchmark)
microbenchmark(HHqr(X),
               CPHHqr(X),
               FPHHqr(X))
*/

sessionInfo () 

R> sessionInfo()
R version 4.0.0 (2020-04-24)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 19.10

Matrix products: default
BLAS:   /usr/lib/x86_64-linux-gnu/openblas/libblas.so.3
LAPACK: /usr/lib/x86_64-linux-gnu/libopenblasp-r0.3.7.so

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] microbenchmark_1.4-7 fBasics_3042.89.1    timeSeries_3062.100 
[4] timeDate_3043.102    Rcpp_1.0.4.11       

loaded via a namespace (and not attached):
 [1] codetools_0.2-16    mvtnorm_1.1-0       lattice_0.20-41    
 [4] zoo_1.8-8           MASS_7.3-51.6       grid_4.0.0         
 [7] spatial_7.3-12      multcomp_1.4-13     Matrix_1.2-18      
[10] fortunes_1.5-4      sandwich_2.5-1      TH.data_1.0-10     
[13] splines_4.0.0       tools_4.0.0         RcppEigen_0.3.3.7.0
[16] survival_3.1-12     parallel_4.0.0      compiler_4.0.0     
[19] colorout_1.2-2     
R> 
R>
...