Как мога да ускоря функцията си с помощта на библиотеката Eigen в C++?

Опитвам се да получа поредица от остатъчна сума от квадрати (RSS) от C++ програма, използвайки for цикъл. И използвах RcppEigen.package.skeleton() за безпроблемно комбиниране на C++ и R. Докато, когато стартирам данни X със 788реда*857кола и Y със 788реда*1кола, времето за изпълнение на C++ програмата е потребителско (4,62 секунди) системно (3,87 секунди) изминало (8,51 s), а времето за изпълнение на програмата R е изминало потребителско (8,68 s) системно (1,78 s) (10,53 s). Програмата C++ е по-малко по-бърза от R. Платформата, която използвах, е win7(X64) с 8G RAM. Как мога да ускоря програмата си? Всяка помощ ще бъде оценена.

Ето програмата на 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);
}