c++ - 如何使用 C++ 中的 Eigen 库加速我的函数？

Question

我正在尝试使用 for 循环从 C++ 程序中获取一系列 Squarts 残差和 (RSS)。我使用 RcppEigen.package.skeleton() 无缝组合 C++ 和 R。而当我运行数据 X 和 788rows*857cols 和 Y 和 788rows*1cols 时，C++ 程序的运行时间是user(4.62s) system(3.87s) elapsed(8.51s)，R程序的运行时间是user(8.68s) system(1.78s) elapsed(10.53s)。 C++ 程序没有 R 快。我使用的平台是 win7(X64) 和 8G 内存。我怎样才能加快我的程序？任何帮助将不胜感激。

这是C++程序:

#include <RcppEigen.h>

//*---get Residual Sum of Squarts via Matrix Operation
//fastLm()
double getRSS(const Eigen::MatrixXd& X, const Eigen::MatrixXd& Y){
   Eigen::MatrixXd RSS=((Y-X*((X.transpose()*X).inverse()*X.transpose()*Y)).transpose())*(Y-X*((X.transpose()*X).inverse()*X.transpose()*Y));
   double RSSd = RSS.determinant();   
   return RSSd;             
}

//*---get F value from RSS and df
double getFval(double RSS1,double RSS2, int n1,int n2,int nObs){
  return (RSS1-RSS2)/(n1-n2)/(RSS2/(nObs-n2-1));      
}

//*---remove p columns from  i-th collumn of matrix
Eigen::MatrixXd removeColumn(const Eigen::MatrixXd& matrix, unsigned int i,int p){
    unsigned int numRows = matrix.rows();
    unsigned int numCols = matrix.cols()-p;

    Eigen::MatrixXd X;
    X=matrix;
    if( i < numCols )
        X.block(0,i,numRows,numCols-i) = matrix.block(0,i+p,numRows,numCols-i);

    X.conservativeResize(numRows,numCols);
    return X;
}

// [[Rcpp::export]]
Rcpp::List getPIcvalue(bool findIn,int p,int n, const Eigen::VectorXd& varIn, const Eigen::MatrixXd& Y,const Eigen::MatrixXd& Xf,const Eigen::MatrixXd& X0){
          //  varIn=(0,1,0,1...,0); p=1 :addition or elimination column; findIn=false,add 1 column of Xf to X0, findIn=false,eliminate 1 column to X0. n=X0.rows();
    bool valid;     
    valid=true;  
    double FitStat1;
    FitStat1 = 1e+10;              

    int pointer;
    pointer=-2;
    double FitStat;
    int nR = n-X0.cols();   // n is the X0.rows()
    int nF;     //nF=nR-1  //findIn=false
    double RSSr;
    double RSSf;
    double F_value;
    RSSr = getRSS(X0,Y);
    int k;
    if(false==findIn){
        k = p;                  
    }else{
        k = -p;      
    }
    Eigen::MatrixXd X(n,X0.cols()+k); 

    if(false==findIn){
        for(int i=0;i<Xf.cols();i++){
            if(0==varIn[i]){
                X<<X0,Xf.col(i);   // X: combine X0 and ith column of Xf                  
                nF = n-X.cols();     
                RSSf = getRSS(X,Y);
                FitStat = getFval(RSSr,RSSf,X.cols(),X0.cols(),n);
                //FitStat = getPvalue(F_value,nF,nR); 
                if(FitStat<FitStat1){
                    FitStat1=FitStat;
                    pointer=i;                    
                }                 
            }//varIn     
        }//for i                 
    }else{
        for(int i=1;i<X0.cols();i++){
            X =  removeColumn(X0,i,p);       
            RSSf = getRSS(X,Y);
            FitStat = getFval(RSSf,RSSr,X0.cols(),X.cols(),n);
            //FitStat = getPvalue(F_value,nR,nF); 
            if(FitStat<FitStat1){
                FitStat1=FitStat;
                pointer=i;                    
            }                 
        }//for i    
    }//findIn 
    return Rcpp::List::create(Rcpp::Named("keyV")=FitStat1,
                              Rcpp::Named("keyP")=pointer+1,
                              Rcpp::Named("keyR")=valid);
}

Answer 1

您对 RSS 矩阵公式的表达效率极低。你这样做:

Eigen::MatrixXd RSS = (
  (Y - X * 
    ( ( X.transpose() * X ).inverse() * X.transpose() * Y ) 
  ).transpose() ) * 
  ( Y - X * 
    ( ( X.transpose() * X ).inverse() * X.transpose() * Y ) 
  );

这显然是非常重复的，并且多次重新计算相同的昂贵操作。转置矩阵应该非常便宜，除非它最终需要一个拷贝。但是反转矩阵(即使是对称正定矩阵，就像这里的情况一样，除非你告诉它，否则 Eigen 无法知道)非常昂贵。哎呀.. 甚至矩阵乘法也很昂贵。

您可能认为 Eigen 做了一些幕后魔术来消除冗余操作并找到最有效的操作序列来获得结果。但是 Eigen 在这方面仍然相当保守(依赖于在编译时解析的保守表达式模板，而实际上它应该使用运行时表达式优化)。所以，它真的不会在这里做那么多。您需要通过自己完成这项工作来帮助它删除冗余操作。

最后，您可以通过执行线性系统解决方案来组合求逆和乘法(而不是 A = inv(X) * B，您可以执行 solve(X * A = B))，它还允许您指定最合适的分解(这里，它是 llt 或 ldlt，具体取决于您期望矩阵 (Xt*X)成为)。

你明白了:

auto Xt = X.transpose(); //<- deduce the type with 'auto' to avoid copy-evaluation of the transpose.
const Eigen::MatrixXd A = X * ( Xt * X ).ldlt().solve(Xt);
const Eigen::MatrixXd Y_AY = Y - A * Y;
Eigen::MatrixXd RSS = Y_AY.transpose() * Y_AY;

但实际上，您可以通过意识到 X * (Xt * X)^-1 * Xt * Y 实际上等同于 X * B 来进一步优化它其中 B 是 X*B = Y 的最小二乘解。如果你使用 QR 方法(不要在这里使用 SVD，它完全矫枉过正而且非常慢，我不明白为什么 Eigen 文档中甚至提到它是线性最小二乘的可行方法(可能是因为 Eigen 人是业余爱好者!))，你可以这样做:

const Eigen::MatrixXd B = X.colPivHouseholderQr().solve( Y );
const Eigen::MatrixXd Y_XB = Y - X * B;
Eigen::MatrixXd RSS = Y_XB.transpose() * Y_XB;

这应该比您以前的速度快得多(至少，就时间复杂度而言，这应该快几个数量级)。

此外，如果 Y 恰好是一个方阵，那么您应该计算 Y_XB 的行列式并将其平方，而不是计算其乘积的行列式自己的转置。这将删除一个矩阵乘法(并复制到 RSS)。

最后，我没有过多研究您的其他函数(调用 getRSS)，但您应该尽一切可能避免重新计算(在每次迭代中)不会改变或不会改变太多的东西，就像 X 的 QR 分解一样。有一些方法可以在 X 的整个变化过程中保持 QR 分解，但这超出了我在这里详细说明的范围，而且可能不是你可以用 Eigen 做的事情。