Как я могу создать подматрицы

Я понимаю, что могу извлекать подматрицы из уже созданной матрицы, но я хочу иметь возможность сначала создавать подматрицы, а затем объединять созданные подматрицы, чтобы сформировать большую матрицу, чтобы сэкономить место и время. Например, в моем примере я хочу иметь возможность создать подматрицу для идентификаторов с NA (1-10) и идентификаторов без NA (11-20), а затем объединить две матрицы вместе, чтобы сформировать большую матрицу, но я не получаю ее , хотел бы, если кто-нибудь может предложить, что должно быть в моих кодах, учитывая, что я буду делать одинаковые расчеты как с NA, так и без NA.

PS: я также хочу иметь возможность сохранять эти подматрицы отдельно, прежде чем объединять их вместе в единую матрицу (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)

r sparse-matrix submatrix
person Viktor    schedule 29.09.2019    source источник
comment
Если вы используете разреженные матрицы, не будет ли проще создать всю матрицу за один раз? Требуются только значения и их позиции в матрице.   -  person Oliver    schedule 29.09.2019
comment
Да, это была моя идея, но для гораздо больших данных (> 400 КБ) он работает уже 2 недели, что довольно медленно, поэтому я спрашиваю об этом и как я могу узнать о значениях и их позициях?   -  person Viktor    schedule 29.09.2019
comment
Это звучит как проблема в коде для генерации разреженной матрицы (как вы это делаете?), больше, чем что-то, что можно было бы исправить с помощью меньших подматриц. Предполагая, что ваши данные имеют 1 миллион переменных, я предполагаю, что матрица, по крайней мере, имеет 1 ненулевое значение в каждой строке и столбце, это не должно занять больше нескольких секунд для создания с использованием пакета Matrix.   -  person Oliver    schedule 30.09.2019
comment
пример: n <- 1e6;d <- rnorm(n);r <- seq(n);c <- sample(r);system.time(mm <- sparseMatrix(i = r, j = c, x = d)). занимает примерно 0,24 секунды на моем маленьком ноутбуке.   -  person Oliver    schedule 30.09.2019
comment
Это казалось знакомым. Я думаю, что использование разреженных матриц в конечном итоге вредит - все ваши NA результаты оказываются неразреженными. См. также: stackoverflow.com/questions/57301390/   -  person Cole    schedule 30.09.2019
comment
@ Оливье, я думаю, почему он медленный из-за ненулевых значений, и он должен вычислять для каждого идентификатора (в соответствии с кодами), поэтому у меня возникла идея вычислить с / без NA, а затем объединить их вместе в конце   -  person Viktor    schedule 30.09.2019
comment
@Cole, это не то же самое, потому что эта связанная страница приведет к сбою R, она не удобна ни по времени, ни по памяти.   -  person Viktor    schedule 30.09.2019
comment
Вы можете изучить bigmemory матрицы (cran.r-project.org /web/packages/bigmemory/index.html), если у вас происходит сбой или заканчивается нехватка памяти при попытке сгенерировать большие матрицы.   -  person Richard J. Acton    schedule 30.09.2019
comment
сейчас речь идет не о нехватке памяти, а о времени/эффективности создания матрицы, поэтому я подумал спросить, как создавать подматрицы, а затем связать их вместе, чтобы в конце сделать большую матрицу.   -  person Viktor    schedule 30.09.2019
comment
Я не думаю, что вы можете сначала сделать подматрицы - последующие результаты зависят от предыдущих результатов в большей матрице. Вы можете просто использовать lower.tri и diag, что уменьшит требования к памяти на 50%. Если вы зададите еще один вопрос с исходным примером из 4000 строк, мне было бы интересно протестировать Rcpp.   -  person Cole    schedule 01.10.2019
comment
Привет, Коул, это ссылка на больший набор данных: ufile.io/irrolcfo спасибо   -  person Viktor    schedule 01.10.2019

Ответы (2)


arrow_upward
0
arrow_downward

Вот решение для 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)
person Cole    schedule 02.10.2019
comment
Пока спасибо, Коля! Я пытаюсь сделать его разреженным, все еще запуская коды на гораздо больших данных: заменил rowSums на rowsum и обновлю вас. - person Viktor; 02.10.2019
comment
@Viktor - см. редактирование. Я перевел ваше исходное решение на Rcpp. На выходе получается разреженная матрица, которая вычисляет ваш больший набор данных за 1 секунду на моем компьютере. - person Cole; 16.10.2019
comment
Как поживаешь? Спасибо!!! Это прекрасно, и я не знаю, смогу ли я добавить ваше имя в мой список благодарностей? - person Viktor; 16.10.2019
comment
@Viktor У меня все отлично. Коул Миллер — мое полное имя, и я был бы очень признателен за подтверждение. - person Cole; 17.10.2019
comment
Я отправлю вам копию моего признания - person Viktor; 17.10.2019

arrow_upward
0
arrow_downward

Чтобы сделать написанную Коулом функцию разреженной, мне пришлось исправить ее с помощью A[t, vec]<- 0.5 * Matrix::rowSums(cbind(A[vec,fam[t,"dad"]],A[vec,fam[t,"mum"]]), na.rm=T)

Спасибо, мы не смогли создать подматрицы, но думаю, что у нас получилось лучше.

person Viktor    schedule 03.10.2019