全透视图QR在C ++中比非透视图QR更快？ （使用Eigen库）

问题描述 投票：0回答：1

通过比较C ++中不同Eigen的QR实现，我得到了怪异的性能结果[我正在通过RcppEigen包从R中访问C ++环境-但是字符串src是只是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.5倍于无旋转QR快两倍？谢谢。

c++ eigen rcpp eigen3
1个回答
0
投票
这是否是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> 

我的代码在下面，我的sessionInfo()也在下面。
代码
#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>
最新问题
© www.soinside.com 2019 - 2024. All rights reserved.