r-我如何使用C++中的特征库来加速我的函数
r - How could I speed up my function using Eigen library in C++?
我正试图使用for循环从C++程序中获得一系列残差平方和(RSS)。我用RcppEigen.package.skeleton()
把C++和R无缝地结合在一起。而当我用788rows*857cols和788rows*1cols运行数据X和Y时,C++程序的运行时间是用户(4.62s)系统(3.87s)运行时间(8.51s),R程序的运行时间是用户(8.68s)系统(1.78s)系统运行时间(10.53s)。C++程序比R快。我使用的平台是带有8G RAM的win7(X64)。我怎样才能加快我的程序?任何帮助都将不胜感激。
这是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);
}
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分解,也许你不能用艾根做点什么。
- 我的神经网络不起作用 [XOR 问题]
- C++我的数学有什么问题,为什么我的代码不能正确循环
- 我的字符计数代码计算错误.为什么
- 为什么我的C#代码在调用回C++COM直到Task时会暂停.等待/线程.加入
- cmake在我的项目中所需的所有静态库都不成功
- 为什么我的代码在输出中增加了93天
- 我的简单if-else语句是如何无法访问的代码
- 为什么我的for循环不能正确获取argv
- 我的项目不会像"undefined reference to `grpc::g_core_codegen_interface'"那样使用未定义的引用错误进行编译
- 0-1背包代码中的错误.我的代码中有什么错误
- 当我的阵列太大时出现分段错误
- 如何确认我的constexpr表达式实际上已经在编译时执行
- 为什么二进制搜索在我的测试中不起作用
- 如何指定我希望我的LIB链接到的DLL文件?-Visual Studio 2019
- 我的代码中有错误吗?使用BGI图形的C++代码对我不起作用
- 当我在main中声明了我的2d数组时,为什么我的程序会退出
- 为什么我的程序在使用预留后没有加速?
- r-我如何使用C++中的特征库来加速我的函数
- 如何加速我的.bmp类
- 如果瓶颈是一个大的矩阵乘法,我可以使用mex文件加速我的MATLAB代码吗?