使用列表以有效的方式在 C++ 中使用 RCPP 返回一堆矩阵
Returning bunch of matrices using RCPP in C++ in an efficient way using a list
我正在尝试 return 使用 RCPP 的一堆矩阵。我下面的代码效率极低。我想知道下面的代码能不能高效
#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]
// [[Rcpp::export]]
Rcpp::List hello(
const arma::rowvec& g,
const int& n,
const int& p,
const arma::mat& S,
const arma::mat& zc,
const arma::rowvec& dl){
Rcpp::List ht(n);
for(int t=0; t < n;++t){
arma::mat hhat(p,n);
hhat.fill(0.0);
for(int i = 0;i < n; ++i){
arma::mat h(p,1);
h.fill(0.0);
if (t > i){
for(int u=i;u <= t; ++u){
arma::rowvec zr = zc.rows(i,i);
h += exp(arma::as_scalar(g*zr.t())) * (zr.t() - S.cols(u,u))*dl(u);
}
}
hhat.cols(i,i) = h;
}
ht[t] = hhat;
}
// Specify list length
Rcpp::List res(1);
res[0] = ht;
return(res);
}
这是例子。
g=c(1,2.1,3.1)
n=1600
p=3
S = matrix(rnorm(4800),nrow=3,ncol=1600)
dl=runif(1600)
z=matrix(runif(4800),nrow=1600,ncol=3)
ptm=proc.time();kkk= hello(g=g,n=n,p=p,S = S,zc=z,dl = dl);proc.time()-ptm;
user system elapsed
31.25 0.00 31.30
如有任何帮助，我们将不胜感激。
遵循更新后的代码。最初我是 returning 一个列表的列表。现在它 return 是一个列表。这将计算时间减少了 10 秒。我希望这段代码可以进一步改进。
#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]
// [[Rcpp::export]]
Rcpp::List hello(
const arma::rowvec& g,
const int& n,
const int& p,
const arma::mat& S,
const arma::mat& zc,
const arma::rowvec& dl){
Rcpp::List ht(n);
for(int t=0; t < n;++t){
arma::mat hhat(p,n);
hhat.zeros();
for(int i = 0;i < n; ++i){
arma::mat h(p,1);
// h.fill(0.0);
h.zeros();
if (t > i){
for(int u=i;u <= t; ++u){
//arma::rowvec zr = zc.rows(i,i);
h += exp(arma::as_scalar(g*zc.row(i).t())) * (zc.row(i).t() - S.col(u))*dl(u);
}
}
hhat.col(i) = h;
}
ht[t] = hhat;
}
// Specify list length
// Rcpp::List res(1);
// res[0] = ht;
return(ht);
}
下面给出了我要实现的公式。
你的问题标题让人认为你看到了将数据返回给 R 的问题。请放心，这不是问题。您可以通过调用一个函数来轻松检查这一点，该函数 returns 所需大小的零矩阵：
#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]
// [[Rcpp::export]]
Rcpp::List minimal(
const arma::rowvec& g,
const int& n,
const int& p,
const arma::mat& S,
const arma::mat& zc,
const arma::rowvec& dl){
Rcpp::List ht(n);
for(int t=0; t < n;++t){
arma::mat hhat(p,n);
hhat.zeros();
ht[t] = hhat;
}
return(ht);
}
在我的系统上，此函数处理您的输入数据大约需要 0.01 秒。换句话说，您的真实函数大部分时间都花在计算实际结果上。
至于优化那部分，如果您能提供您正在尝试实现的内容的想法将会很有帮助，例如借助数学公式。目前，我只能做一些简单的更改：
在 i 循环中，您只为 t > i 做一些事情。因此，让循环 运行 直到 i < t.
u 循环可以表述为一个 matrix-vector 产品，对此存在有效的实现。
有了这样的改变，我最终得到了
#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]
// [[Rcpp::export]]
Rcpp::List hello(
const arma::rowvec& g,
const int& n,
const int& p,
const arma::mat& S,
const arma::mat& zc,
const arma::rowvec& dl){
Rcpp::List ht(n);
for(int t=0; t < n;++t){
arma::mat hhat(p,n);
hhat.zeros();
for(int i = 0;i < t; ++i){
arma::mat Sit = S.cols(i,t);
hhat.col(i) = - exp(arma::as_scalar(g*zc.row(i).t())) *
(Sit.each_col() - zc.row(i).t()) * dl.subvec(i,t).t();
}
ht[t] = hhat;
}
return(ht);
}
在我的系统上，这比您的代码快两倍。不过，它可能会变得更快。
在我的另一个回答中，我查看了 returning 数据的效率和简单的优化。在这里我想看一些不同的东西：算法的优化。
您想为 0 <= i, t < n 和 i < t 计算 hhat(i, t)。查看您的公式，我们发现 hhat 对 i 和 t 的依赖性非常不同。特别地，hhat(i, t + 1) 可以写成 hhat(i, t) + something。现在你的外循环已经结束 t 并且你正在重新计算所有这些中间值。通过切换循环顺序，很容易将每个此类计算仅执行一次，从而将算法简化为两个嵌套循环。这意味着您必须单独生成结果矩阵。因为你不能在 Rcpp::List 中存储 arma::mat，我需要一个额外的 std::vector 来存储：
// [[Rcpp::depends(RcppArmadillo)]]
#include <RcppArmadillo.h>
// [[Rcpp::export]]
Rcpp::List hello(
const arma::rowvec& g,
const int& n,
const int& p,
const arma::mat& S,
const arma::mat& zc,
const arma::rowvec& dl){
std::vector<arma::mat> foo(n);
for(int t=0; t < n;++t){
arma::mat hhat(p,n);
hhat.zeros();
foo[t] = hhat;
}
for(int i = 0;i < n; ++i){
arma::mat h = exp(arma::as_scalar(g*zc.row(i).t())) * (zc.row(i).t() - S.col(i))*dl(i);
for(int t=i+1; t < n;++t){
h += exp(arma::as_scalar(g*zc.row(i).t())) * (zc.row(i).t() - S.col(t))*dl(t);
foo[t].col(i) = h;
}
}
Rcpp::List ht(n);
for(int t=0; t < n;++t){
ht[t] = foo[t];
}
return(ht);
}
// [[Rcpp::export]]
Rcpp::List hello_orig(
const arma::rowvec& g,
const int& n,
const int& p,
const arma::mat& S,
const arma::mat& zc,
const arma::rowvec& dl){
Rcpp::List ht(n);
for(int t=0; t < n;++t){
arma::mat hhat(p,n);
hhat.zeros();
for(int i = 0;i < n; ++i){
arma::mat h(p,1);
h.zeros();
if (t > i){
for(int u=i;u <= t; ++u){
h += exp(arma::as_scalar(g*zc.row(i).t())) * (zc.row(i).t() - S.col(u))*dl(u);
}
}
hhat.col(i) = h;
}
ht[t] = hhat;
}
return(ht);
}
/***R
g=c(1,2.1,3.1)
n=1600
p=3
S = matrix(rnorm(p*n),nrow=p,ncol=n)
dl=runif(n)
z=matrix(runif(p*n),nrow=n,ncol=p)
bench::mark(hello_orig(g=g,n=n,p=p,S = S,zc=z,dl = dl),
hello(g=g,n=n,p=p,S = S,zc=z,dl = dl))
*/
结果：
# A tibble: 2 x 13
expression min median `itr/sec` mem_alloc
<bch:expr> <bch:> <bch:> <dbl> <bch:byt>
1 hello_orig(g = g, n = n, p = p, S = S, zc = z, dl = dl) 14.2s 14.2s 0.0703 58.7MB
2 hello(g = g, n = n, p = p, S = S, zc = z, dl = dl) 53.9ms 85.9ms 11.1 58.7MB
# … with 8 more variables: `gc/sec` <dbl>, n_itr <int>, n_gc <dbl>, total_time <bch:tm>,
# result <list>, memory <list>, time <list>, gc <list>
速度提高 100 倍以上！
您可以通过在评论中使用 @coatless 的建议使用 arma::cube 来获得更清晰（甚至可能更快的代码）。不过，最紧凑的形式将为您提供不同的 return 结构。您将得到一个 p x n x n 数组，而不是 p x n 列表：
// [[Rcpp::depends(RcppArmadillo)]]
#include <RcppArmadillo.h>
// [[Rcpp::export]]
arma::cube coatless(
const arma::rowvec& g,
const int& n,
const int& p,
const arma::mat& S,
const arma::mat& zc,
const arma::rowvec& dl){
arma::cube ht(p, n, n);
ht.zeros();
for(int i = 0;i < n; ++i){
arma::mat h = exp(arma::as_scalar(g*zc.row(i).t())) * (zc.row(i).t() - S.col(i))*dl(i);
for(int t=i+1; t < n;++t){
h += exp(arma::as_scalar(g*zc.row(i).t())) * (zc.row(i).t() - S.col(t))*dl(t);
ht.slice(t).col(i) = h;
}
}
return(ht);
}
我正在尝试 return 使用 RCPP 的一堆矩阵。我下面的代码效率极低。我想知道下面的代码能不能高效
#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]
// [[Rcpp::export]]
Rcpp::List hello(
const arma::rowvec& g,
const int& n,
const int& p,
const arma::mat& S,
const arma::mat& zc,
const arma::rowvec& dl){
Rcpp::List ht(n);
for(int t=0; t < n;++t){
arma::mat hhat(p,n);
hhat.fill(0.0);
for(int i = 0;i < n; ++i){
arma::mat h(p,1);
h.fill(0.0);
if (t > i){
for(int u=i;u <= t; ++u){
arma::rowvec zr = zc.rows(i,i);
h += exp(arma::as_scalar(g*zr.t())) * (zr.t() - S.cols(u,u))*dl(u);
}
}
hhat.cols(i,i) = h;
}
ht[t] = hhat;
}
// Specify list length
Rcpp::List res(1);
res[0] = ht;
return(res);
}
这是例子。
g=c(1,2.1,3.1)
n=1600
p=3
S = matrix(rnorm(4800),nrow=3,ncol=1600)
dl=runif(1600)
z=matrix(runif(4800),nrow=1600,ncol=3)
ptm=proc.time();kkk= hello(g=g,n=n,p=p,S = S,zc=z,dl = dl);proc.time()-ptm;
user system elapsed
31.25 0.00 31.30
如有任何帮助，我们将不胜感激。
遵循更新后的代码。最初我是 returning 一个列表的列表。现在它 return 是一个列表。这将计算时间减少了 10 秒。我希望这段代码可以进一步改进。
#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]
// [[Rcpp::export]]
Rcpp::List hello(
const arma::rowvec& g,
const int& n,
const int& p,
const arma::mat& S,
const arma::mat& zc,
const arma::rowvec& dl){
Rcpp::List ht(n);
for(int t=0; t < n;++t){
arma::mat hhat(p,n);
hhat.zeros();
for(int i = 0;i < n; ++i){
arma::mat h(p,1);
// h.fill(0.0);
h.zeros();
if (t > i){
for(int u=i;u <= t; ++u){
//arma::rowvec zr = zc.rows(i,i);
h += exp(arma::as_scalar(g*zc.row(i).t())) * (zc.row(i).t() - S.col(u))*dl(u);
}
}
hhat.col(i) = h;
}
ht[t] = hhat;
}
// Specify list length
// Rcpp::List res(1);
// res[0] = ht;
return(ht);
}
下面给出了我要实现的公式。
你的问题标题让人认为你看到了将数据返回给 R 的问题。请放心，这不是问题。您可以通过调用一个函数来轻松检查这一点，该函数 returns 所需大小的零矩阵：
#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]
// [[Rcpp::export]]
Rcpp::List minimal(
const arma::rowvec& g,
const int& n,
const int& p,
const arma::mat& S,
const arma::mat& zc,
const arma::rowvec& dl){
Rcpp::List ht(n);
for(int t=0; t < n;++t){
arma::mat hhat(p,n);
hhat.zeros();
ht[t] = hhat;
}
return(ht);
}
在我的系统上，此函数处理您的输入数据大约需要 0.01 秒。换句话说，您的真实函数大部分时间都花在计算实际结果上。
至于优化那部分，如果您能提供您正在尝试实现的内容的想法将会很有帮助，例如借助数学公式。目前，我只能做一些简单的更改：
在
i循环中，您只为t > i做一些事情。因此，让循环 运行 直到i < t. 就足够了
u循环可以表述为一个 matrix-vector 产品，对此存在有效的实现。
有了这样的改变，我最终得到了
#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]
// [[Rcpp::export]]
Rcpp::List hello(
const arma::rowvec& g,
const int& n,
const int& p,
const arma::mat& S,
const arma::mat& zc,
const arma::rowvec& dl){
Rcpp::List ht(n);
for(int t=0; t < n;++t){
arma::mat hhat(p,n);
hhat.zeros();
for(int i = 0;i < t; ++i){
arma::mat Sit = S.cols(i,t);
hhat.col(i) = - exp(arma::as_scalar(g*zc.row(i).t())) *
(Sit.each_col() - zc.row(i).t()) * dl.subvec(i,t).t();
}
ht[t] = hhat;
}
return(ht);
}
在我的系统上，这比您的代码快两倍。不过，它可能会变得更快。
在我的另一个回答中，我查看了 returning 数据的效率和简单的优化。在这里我想看一些不同的东西：算法的优化。
您想为 0 <= i, t < n 和 i < t 计算 hhat(i, t)。查看您的公式，我们发现 hhat 对 i 和 t 的依赖性非常不同。特别地，hhat(i, t + 1) 可以写成 hhat(i, t) + something。现在你的外循环已经结束 t 并且你正在重新计算所有这些中间值。通过切换循环顺序，很容易将每个此类计算仅执行一次，从而将算法简化为两个嵌套循环。这意味着您必须单独生成结果矩阵。因为你不能在 Rcpp::List 中存储 arma::mat，我需要一个额外的 std::vector 来存储：
// [[Rcpp::depends(RcppArmadillo)]]
#include <RcppArmadillo.h>
// [[Rcpp::export]]
Rcpp::List hello(
const arma::rowvec& g,
const int& n,
const int& p,
const arma::mat& S,
const arma::mat& zc,
const arma::rowvec& dl){
std::vector<arma::mat> foo(n);
for(int t=0; t < n;++t){
arma::mat hhat(p,n);
hhat.zeros();
foo[t] = hhat;
}
for(int i = 0;i < n; ++i){
arma::mat h = exp(arma::as_scalar(g*zc.row(i).t())) * (zc.row(i).t() - S.col(i))*dl(i);
for(int t=i+1; t < n;++t){
h += exp(arma::as_scalar(g*zc.row(i).t())) * (zc.row(i).t() - S.col(t))*dl(t);
foo[t].col(i) = h;
}
}
Rcpp::List ht(n);
for(int t=0; t < n;++t){
ht[t] = foo[t];
}
return(ht);
}
// [[Rcpp::export]]
Rcpp::List hello_orig(
const arma::rowvec& g,
const int& n,
const int& p,
const arma::mat& S,
const arma::mat& zc,
const arma::rowvec& dl){
Rcpp::List ht(n);
for(int t=0; t < n;++t){
arma::mat hhat(p,n);
hhat.zeros();
for(int i = 0;i < n; ++i){
arma::mat h(p,1);
h.zeros();
if (t > i){
for(int u=i;u <= t; ++u){
h += exp(arma::as_scalar(g*zc.row(i).t())) * (zc.row(i).t() - S.col(u))*dl(u);
}
}
hhat.col(i) = h;
}
ht[t] = hhat;
}
return(ht);
}
/***R
g=c(1,2.1,3.1)
n=1600
p=3
S = matrix(rnorm(p*n),nrow=p,ncol=n)
dl=runif(n)
z=matrix(runif(p*n),nrow=n,ncol=p)
bench::mark(hello_orig(g=g,n=n,p=p,S = S,zc=z,dl = dl),
hello(g=g,n=n,p=p,S = S,zc=z,dl = dl))
*/
结果：
# A tibble: 2 x 13
expression min median `itr/sec` mem_alloc
<bch:expr> <bch:> <bch:> <dbl> <bch:byt>
1 hello_orig(g = g, n = n, p = p, S = S, zc = z, dl = dl) 14.2s 14.2s 0.0703 58.7MB
2 hello(g = g, n = n, p = p, S = S, zc = z, dl = dl) 53.9ms 85.9ms 11.1 58.7MB
# … with 8 more variables: `gc/sec` <dbl>, n_itr <int>, n_gc <dbl>, total_time <bch:tm>,
# result <list>, memory <list>, time <list>, gc <list>
速度提高 100 倍以上！
您可以通过在评论中使用 @coatless 的建议使用 arma::cube 来获得更清晰（甚至可能更快的代码）。不过，最紧凑的形式将为您提供不同的 return 结构。您将得到一个 p x n x n 数组，而不是 p x n 列表：
// [[Rcpp::depends(RcppArmadillo)]]
#include <RcppArmadillo.h>
// [[Rcpp::export]]
arma::cube coatless(
const arma::rowvec& g,
const int& n,
const int& p,
const arma::mat& S,
const arma::mat& zc,
const arma::rowvec& dl){
arma::cube ht(p, n, n);
ht.zeros();
for(int i = 0;i < n; ++i){
arma::mat h = exp(arma::as_scalar(g*zc.row(i).t())) * (zc.row(i).t() - S.col(i))*dl(i);
for(int t=i+1; t < n;++t){
h += exp(arma::as_scalar(g*zc.row(i).t())) * (zc.row(i).t() - S.col(t))*dl(t);
ht.slice(t).col(i) = h;
}
}
return(ht);
}