Более эффективная функция для генерации специальной матрицы в R

Question

Я пытаюсь запрограммировать функцию R для создания специальной матрицы с именем Ui в этой статье (стр. 5): http://joshuachan.org/papers/Chan-Jeliazkov-2009.pdf

До сих пор я разработал эту функцию, котораяЯ подозреваю, что может быть улучшено (то есть сделано более эффективным):

create_Uj <- function(uj) {
  q <- length(uj) 
  if (q == 1) return(0)
  nr <- q
  nc <- q*(q-1)/2
  Uj <- matrix(0, nrow = nr, ncol = nc)
  for (kk in 2:nr) {
    uj_sub <- uj[1:(kk-1)]
    Uj[kk, 1:(kk*(kk-1)/2)] <- c(rep(0, (kk-1)*(kk-2)/2), uj[1:(kk-1)])
  }
  -Uj
}

create_Uj(1:4)

     [,1] [,2] [,3] [,4] [,5] [,6]
[1,]    0    0    0    0    0    0
[2,]   -1    0    0    0    0    0
[3,]    0   -1   -2    0    0    0
[4,]    0    0    0   -1   -2   -3

Есть ли лучший способ закодировать это?

Примечание. Я использую индекс j вместо i в статье

Answer 1

Можно также написать функцию, используя вложенные циклы:

create_Uj3 = function(uj){
  nr <- length(uj) 
  if (nr == 1){
    return(0)
  } 
  nc <- nr*(nr-1)/2
  Uj <- matrix(0, nrow = nr, ncol = nc)

  for (kk in 2:nr) {
    for (ll in 1:(kk-1)){
      Uj[kk, ((kk-1)*(kk-2)/2) + ll] <- uj[ll]
    }
  }
  return(-Uj)
}

Эквивалент Rcpp:

library(Rcpp)
cppFunction('NumericMatrix create_Uj_rcpp(NumericVector uj) {
  const int nr = uj.size();
  if(nr == 1){
    return 0;
  }
  const int nc = (nr*(nr-1))/2;
  NumericMatrix Uj = NumericMatrix(nr, nc);

  for(int i = 1; i < nr; i++) {
    for(int j = 0; j <= (i - 1); j++){
      Uj(i,(i*(i-1)/2) + j) = -uj[j];
    }
  }
  return Uj;
}')

Тесты:

> identical(create_Uj(1:300), create_Uj2(1:300))
[1] TRUE
> identical(create_Uj(1:300), create_Uj3(1:300))
[1] TRUE
> identical(create_Uj(1:300), create_Uj_rcpp(1:300))
[1] TRUE

library(microbenchmark)
microbenchmark(create_Uj(1:300), 
               create_Uj2(1:300), 
               create_Uj3(1:300),
               create_Uj_rcpp(1:300), 
               unit = 'relative')

Unit: relative
                  expr      min       lq     mean   median       uq      max neval
      create_Uj(1:300) 6.299113 6.874493 4.351439 5.128115 4.059644 2.105598   100
     create_Uj2(1:300) 2.859025 3.827864 2.524505 2.873346 2.233600 1.618327   100
     create_Uj3(1:300) 3.078410 4.111635 2.552434 3.109537 2.259824 1.418917   100
 create_Uj_rcpp(1:300) 1.000000 1.000000 1.000000 1.000000 1.000000 1.000000   100

create_Uj_rcpp выигрывает по скорости.Обратите внимание, что метод Base R, вложенный в цикл (create_Uj3), немного медленнее, чем решение Райана (create_Uj2), но все же намного быстрее, чем функция OP (create_Uj).

Answer 2

Вы можете получить приличное ускорение, просто избегая создания rep(0, (kk-1)*(kk-2)/2).Тот факт, что удаление этого шага значительно ускоряет процесс, заставляет меня думать, что вряд ли вы станете намного быстрее без использования Rcpp

create_Uj2 <- function(uj) {
  q <- length(uj) 
  if (q == 1) return(0)
  nr <- q
  nc <- q*(q-1)/2
  Uj <- matrix(0, nrow = nr, ncol = nc)
  for (kk in 2:nr) {
    Uj[kk, ((kk-1)*(kk-2)/2 + 1):(kk*(kk-1)/2)] <- uj[1:(kk-1)]
  }
  -Uj
}

all.equal(create_Uj2(1:400), create_Uj(1:400))
# [1] TRUE
microbenchmark(
  create_Uj2(1:400),
  create_Uj(1:400),
  times = 10,
  unit = 'relative'
)
# Unit: relative
#               expr      min       lq     mean   median       uq       max neval
#  create_Uj2(1:400) 1.000000 1.000000 1.000000 1.000000 1.000000 1.0000000    10
#   create_Uj(1:400) 2.070708 1.847242 1.529658 2.028489 1.496275 0.9441935    10