如何创建子矩阵答案

【问题标题】：How can I create submatrices如何创建子矩阵
【发布时间】：2020-01-29 04:02:24
【问题描述】：

我知道我可以从已经创建的矩阵中提取子矩阵，但我希望能够先创建子矩阵，然后组合创建的子矩阵以形成更大的矩阵以节省空间和时间。例如，在我的示例中，我希望能够为具有 NA (1-10) 的 ID 和不具有 NA (11-20) 的 ID 创建一个子矩阵，然后将这两个矩阵组合在一起以形成一个更大的矩阵，但我没有得到它，如果有人可以建议我的代码中应该包含的内容，那么我将对使用 NA 和不使用 NA 进行相同的计算。

P.S：我还希望能够在将这些子矩阵合并到一个奇异矩阵（20x20）之前单独保存它们

dorm<-function(data)
{ 
  library(Matrix)
  n<-max(as.numeric(fam[,"ID"])) 
  t<-min(as.numeric(fam[,"ID"])) 
  A <- sparseMatrix(i = n, j=n, x=n)
  while(t <=n) {

    for( t in t:n ){

      s <- max(fam[t,"dad"],fam[t,"mum"]) 
      d <- min(fam[t,"dad"],fam[t,"mum"])

      if( !is.na(s) ){ 
        if( !is.na(d) ){
          A[t,t] = 2-0.5^(fam[t,"GEN"]-1)+0.5^(fam[t,"GEN"])*A[fam[t,"dad"],fam[t,"mum"]]
          tmp = 0.5 * (A[1:(t-1),s] + A[1:(t-1),d])
          A[t, 1:(t-1)] = tmp
          A[1:(t-1), t] = tmp
        } else {
          A[t,t] = 2-0.5^(fam[t,"GEN"]-1)
          tmp = 0.5 * A[1:(t-1),s]
          A[t, 1:(t-1)] = tmp
          A[1:(t-1), t] = tmp
        }
      } else {
        A[t,t] = 2-0.5^(fam[t,"GEN"]-1)
      }
      message(" MatbyGEN: ", t)
    }

    return(A)
  }
}

fam <- structure(list(ID = c(1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 
11L, 12L, 13L, 14L, 18L, 15L, 16L, 17L, 20L, 19L), dad = c(NA, 
NA, NA, NA, NA, NA, NA, NA, NA, NA, 1L, 1L, 4L, 6L, 4L, 10L, 
12L, 13L, 13L, 14L), mum = c(NA, NA, NA, NA, NA, NA, NA, NA, 
NA, NA, 2L, 3L, 2L, 5L, 11L, 11L, 5L, 3L, 7L, 2L), GEN = c(1L, 
1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 
3L, 3L, 3L)), class = "data.frame", row.names = c(NA, -20L))

A <- dorm(fam)

【问题讨论】：

如果你使用稀疏矩阵，一次性创建整个矩阵不是更简单吗？它只需要矩阵中的值及其位置。
是的，这是我的想法，但对于更大的数据 (>400k)，它已经运行了 2 周，这非常慢，这就是我要求这个的原因以及我该如何解决价值观和他们的立场？
这听起来像是生成稀疏矩阵的代码中的一个问题（你是怎么做到的？），而不是使用较小的子矩阵来解决的问题。假设您的数据有 1M 变量，我假设矩阵在每一行和每一列中至少有 1 个非零值，使用 Matrix 包创建这应该不会超过几秒钟。
例如：n <- 1e6;d <- rnorm(n);r <- seq(n);c <- sample(r);system.time(mm <- sparseMatrix(i = r, j = c, x = d))。在我的小笔记本电脑上大约需要 0.24 秒。
这似乎很熟悉。我认为使用稀疏矩阵最终会受到伤害 - 你所有的 NA 结果结果都是非稀疏的。另见：stackoverflow.com/questions/57301390/…

标签： r sparse-matrix submatrix

【解决方案1】：

这是一个rcpp 解决方案。它在大型数据集上快了约 50 倍（1 秒对 50 秒）：

#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]
using namespace Rcpp;
using namespace arma;

// [[Rcpp::export]]
sp_mat rcpp_dorm_sp(IntegerVector ID, IntegerVector dad, IntegerVector mum, IntegerVector gen){
  int n; 
  int s; int d;

  double tmp;

  sp_mat A(dad.size(), dad.size());

  A.diag().ones();
  n = max(ID); 

  for(int t = 0; t < n; t++){
    s = std::max(dad[t], mum[t]); 
    d = std::min(dad[t], mum[t]);

    A(t,t) = 2-pow(0.5, gen[t] - 1);

    if ((s>0) & (d>0) ) { 
      A(t,t) +=  pow(0.5, gen[t])*A(dad[t]-1,mum[t]-1);
      for(int j = 0; j < t; j++){

        tmp = 0.5 * (A(j, dad[t]-1) + A(j, mum[t]-1));
        if (tmp > 0){
          A(t,j) = tmp;
          A(j,t) = tmp;
        }
      }
    } else if ((s>0) & (d==0)) {

      for(int j = 0; j < t; j++){
        tmp = 0.5 * A(j, s-1);
        if (tmp > 0){
          A(t,j) = tmp;
          A(j,t) = tmp;
        }
      }
    }
  }

  return(A);
}

还有R 部分：

fam_mid <- structure(list(ID = c(1L, 2L, 3L, 4L, 5L, 6L, 7L, 8L, 9L, 10L, 
                                         11L, 12L, 13L, 14L, 18L, 15L, 16L, 17L, 20L, 19L),
                                  dad = c(NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 1L, 1L, 4L, 6L, 4L, 10L, 
                                          12L, 13L, 13L, 14L),
                                  mum = c(NA, NA, NA, NA, NA, NA, NA, NA, NA, NA, 2L, 3L, 2L, 5L, 11L, 11L, 5L, 3L, 7L, 2L)
                                  , GEN = c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 
                                            3L, 3L, 3L)), class = "data.frame", row.names = c(NA, -20L))

rcpp_dorm_sp(fam_cpp$ID, fam_cpp$dad, fam_cpp$mum, fam_cpp$GEN)

【讨论】：

到目前为止，谢谢，科尔！我正在尝试使其稀疏，仍然在更大的数据上运行代码：将 rowSums 替换为 rowsum 并将更新您
@Viktor - 见编辑。我将您的原始解决方案翻译成Rcpp。输出是一个稀疏矩阵，并在我的计算机上在 1 秒内计算出你更大的数据集。
你好吗？谢谢！！！这很漂亮，我不知道我是否可以将你的名字添加到我的致谢列表中？
@Viktor 我做得很好。 Cole Miller 是我的全名，非常感谢您的认可
我会将我的确认书副本发送给您

【解决方案2】：

为了使 Cole 编写的函数变得稀疏，我不得不使用 A[t, vec]<- 0.5 * Matrix::rowSums(cbind(A[vec,fam[t,"dad"]],A[vec,fam[t,"mum"]]), na.rm=T) 修复它

到目前为止，谢谢，我们无法创建子矩阵，但认为我们做得更好

【讨论】：