矩阵列的操作

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

送您一个全额奖学金名额~ !

立即领取

感谢您参与论坛问题回答

经管之家送您两个论坛币！

+2 论坛币

如图所示，有没有可以快速得到 39,66,99 的方法，因为当我列很多的时候，相当于做 p*(p-1)/2次相加。

扫码加我 拉你入群

请注明：姓名-公司-职位

以便审核进群资格，未注明则拒绝

分享0 收藏2 回帖

关键词：如图所示 有没有 如图所示

1. x <- matrix(1:9, 3)
   
2. n <- ncol(x)
   
3. colprod <- vector(length = n)
   
4. for (i in 1:(n-1)) {
   
5.   for (j in (i+1):n)
   
6.     colprod <- colprod + x[,i]*x[,j]   
   
7. }
   
8. colprod

复制代码

楼上同学的解决方案中有一个小错误，代码的第三行应该是：
colprod <- vector(length = nrow(x))
可能是编的快，疏忽了。

楼主的这个问题，我尝试用向量化来解决，发现力不从心。该问题似乎是循环起来很简单，但是向量化非常吃力的那种问题，并且，该问题的计算量随着列数的增长而呈现几何级数地增加的。在我的电脑中，当列数过1000后，耗时就比较明显了。于是，我想用Rcpp包的C++程序来解决这个问题：

1. # 需安装和配置Rtools工具
   
2. library(Rcpp)
   
3. cppFunction('NumericVector colprodC(NumericMatrix x) {
   
4.   int nrow = x.nrow(), ncol = x.ncol();
   
5.   NumericVector out(nrow);
   
6. 
   
7.   for (int i = 0; i < (ncol - 1); i++) {
   
8.     for (int j = (i + 1); j < ncol; j++) {
   
9.       for (int r = 0; r < nrow; r++) {
   
10.         out[r] += x(r, i) * x(r, j);
    
11.       }
    
12.     }
    
13.   }
    
14.   return out;
    
15. }')
    
16. 
    
17. colprodC(mat)

复制代码

在我的电脑上运行，速度提升大概36倍左右，不是很让人满意，是不是我的程序还有问题，请大神指点。

stzhao 发表于 2016-3-15 21:28

谢谢。