вычислить только верхний треугольник матрицы

Question

У меня есть вектор:

v1 = c(1,2,3) 

Из этого вектора Я хочу, чтобы создать матрицу, где элемент на i,j положение будет sum членов вектора на i,j позиции:

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

Вопросы:

1. i,j и j,i - это то же самое, поэтому нет смысла вычислять его 2x для лучшей производительности. Как достичь этого?
2. Как создать также вариант, который не будет вычислять элементы, если i == j и просто возвращает NA? Я не прошу команду diag(m) <- NA, я спрашиваю, как предотвратить вычисление этих элементов.

PS: Это сокращенная версия моей проблемы

Answer 1

Существует подход, который гораздо быстрее, чем прямое вычисление с 2 вложенных циклов. Это не оптимизировано в терминах, которые вы описали в вопросе 1, но это довольно быстро, потому что оно векторизовано. Возможно, этого будет достаточно для вашей цели. (или даже матрица) подход Векторизованного само:

f1 <- function(x){ 
    n <- length(x) 
    m <- matrix(rep(x,n),n) 
    m + t(m) 
} 
> f1(1:3) 
     [,1] [,2] [,3] 
[1,] 2 3 4 
[2,] 3 4 5 
[3,] 4 5 6 

Мы также можем создать функцию для простого подхода для выполнения теста. Эта функция делает даже меньше, чем необходимо: она вычисляет только верхний треугольник, но мы увидим, что она намного медленнее.

f2 <- function(x){ 
    n <- length(x) 
    m <- matrix(rep(NA,n^2),n) 
    for(i in 1:(n-1)){ 
     for(j in (i+1):n) m[i,j] <- x[[i]] + x[[j]] 
    } 
    m 
} 
> f2(1:3) 
     [,1] [,2] [,3] 
[1,] NA 3 4 
[2,] NA NA 5 
[3,] NA NA NA 

Benchmark:

library(microbenchmark)  
> microbenchmark(f1(1:100), f2(1:100)) 
Unit: microseconds 
     expr  min   lq  mean median  uq  max neval 
f1(1:100) 124.775 138.6175 181.6401 187.731 196.454 294.301 100 
f2(1:100) 10227.337 10465.1285 11000.1493 10616.830 10907.148 15826.259 100 

Answer 2

Я думаю, что это хороший пример, чтобы узнать некоторые (основные) материал вызова скомпилированного кода внутри R для ускорения. Здесь я предоставлю функцию C syadd (симметричная добавка) и ее R-оберточная функция, называемая также syadd.

/* symmetric add */ 
/* save this file as "syadd.c"*/ 

#include <R.h> 
#include <Rinternals.h> 

SEXP syadd (SEXP N, SEXP V, SEXP compute_diag) { 
    int n = asInteger(N), diag = asInteger(compute_diag); 
    /* initialize matrix X as a NA/NaN matrix */ 
    SEXP X = PROTECT(allocVector(REALSXP, n * n)); 
    double *x = REAL(X), *ptr_x = x, *v_end = x + n * n; 
    while (ptr_x < v_end) *(ptr_x++) = NA_REAL; 
    /* C interface */ 
    double *v = REAL(V), *vj = v, *vi, *vi_end = v, tmp; 
    ptr_x = x; v_end = v + n; 
    if (diag == 1) { 
    /* compute upper triangular (including diagonal) */ 
    while (vj < v_end) { 
     tmp = *vj; 
     ptr_x = x; vi = v; 
     while (vi <= vi_end) { 
     *ptr_x = (*vi) * tmp; 
     vi++; ptr_x++; 
     } 
     x += n; vi_end++; vj++; 
     } 
    } else { 
    /* compute upper triangular only */ 
    while (vj < v_end) { 
     tmp = *vj; 
     ptr_x = x; vi = v; 
     while (vi < vi_end) { 
     *ptr_x = (*vi) * tmp; 
     vi++; ptr_x++; 
     } 
     x += n; vi_end++; vj++; 
     } 
    } 
UNPROTECT(1); 
return X; 
} 

Чтобы использовать этот код в R, вы должны скомпилировать его в общую библиотеку (.so в Linux или .dll в окнах). Я не знаю, как работать вокруг на Windows, но это, как сделать на Linux:

R CMD SHLIB -c syadd.c 

Тогда давайте проверим его в R:

## R wrapper function 
syadd <- function (v, diag = TRUE) { 
    n <- length(v); v <- as.numeric(v); diag <- as.integer(diag) 
    ## load shared library 
    dyn.load("syadd.so") 
    X <- .Call("syadd", n, v, diag) 
    ## add "dim" attribute to get a matrix rather than a vector 
    attr(X, "dim") <- c(n, n) 
    return(X) 
    } 

Если мы устанавливаем diag = FALSE, диагональные элементы не вычисляется. Давайте бенчмарк его от R-х base:::outer:

## benchmarking 
v <- 1:10000 
system.time(syadd(v,diag = TRUE)) 
gc() 
system.time(outer(v,v,"+")) 
gc() 

На моем ноутбуке, он принимает 0.97s для syadd, но 1.91s для outer. Я думаю, вам нужно разместить мусорную сборку gc() здесь, так как иначе вы можете поменять местами, что приводит к замедлению.

Комментарии:

1. Я только показывает основной способ вызова скомпилированного кода, явным образом с помощью dyn.load(). Можно использовать inline пакет cfunction. См. Примеры R's C interface - advanced R.
2. Было бы хорошо, если бы кто-то мог добавить детали реализации в Windows.