Pytanie Efektywne obliczanie macierzy var-covar w R


Szukam przyrostów wydajności w obliczaniu macierzy kowariancji (auto) z poszczególnych pomiarów w czasie t z t, t-1itd ..

W macierzy danych każdy wiersz reprezentuje osobnika, a każda kolumna reprezentuje miesięczne pomiary (kolumny są w kolejności czasu). Podobne do poniższych danych (chociaż z nieco większą ko-wariancją).

# simulate data
set.seed(1)
periods <- 70L
ind <- 90000L
mat <- sapply(rep(ind, periods), rnorm)

Poniżej znajduje się (brzydki) kod, który wymyśliłem, aby uzyskać macierz kowariancji dla pomiarów / opóźnionych pomiarów. Uruchomienie trwa prawie 4 sekundy. Jestem pewien, że przechodząc do data.table, myśląc więcej i nie polegając na pętlach, mógłbym skrócić czas o dużą ilość. Ale ponieważ matryce kowariancji są wszechobecne, podejrzewam, że istnieje już standardowy (i skuteczny) sposób robienia tego w R, o którym powinienem wiedzieć.

# Get variance covariance matrix for 0-5 lags    
n_lags <- 5L # Number of lags
vcov <- matrix(0, nrow = n_lags + 1L, ncol = n_lags + 1)
for (i in 0L:n_lags) {
  for (j in i:n_lags) {
    vcov[j + 1L, i + 1L] <- 
      sum(mat[, (1L + (j - i)):(periods - i)] *
          mat[, 1L:(periods - j)]) /
      (ind * (periods - j) - 1)
  }
}
round(vcov, 3)

       [,1]   [,2]  [,3]  [,4]  [,5]  [,6]
[1,]  1.001  0.000 0.000 0.000 0.000 0.000
[2,]  0.000  1.001 0.000 0.000 0.000 0.000
[3,]  0.000  0.000 1.001 0.000 0.000 0.000
[4,]  0.000  0.000 0.000 1.001 0.000 0.000
[5,] -0.001  0.000 0.000 0.000 1.001 0.000
[6,]  0.000 -0.001 0.000 0.000 0.000 1.001

14
2017-07-11 22:15


pochodzenie


Spójrz na cov() - Andrew Gustar
Dzięki. Ale jeśli sugerujesz cov(mat)[1:6, 1:6] to trochę inaczej ... ponieważ nie szukam kowariancji t=1 z t=2 ale generalnie t z t-1... ale może będę mógł użyć tej funkcji, jeśli skonfiguruję swoją macierz inaczej (?). - snoram
Sprawdź ?ccf funkcjonować? - Ben Bolker
Być może covmat <- cov(mat) Aby uzyskać kowariancje między kolumnami, można wyodrębnić dane liczbowe dotyczące opóźnienia Lz czymś podobnym diag(covmat[-(1:L),]) - Andrew Gustar
Jeśli naprawdę zależy Ci na wydajności, myślę, że twoim najlepszym ujęciem byłby algorytm Rcpp z jednym przebiegiem. Jedyny problem, jaki widzę z tym kodem R, polega na tym, że tworzy on wiele macierzy (kopii). - F. Privé


Odpowiedzi:


@FA. Privé's Rcpp wdrożenie jest dobrym początkiem, ale możemy zrobić to lepiej. W głównym algorytmie dostarczonym przez OP zauważysz, że istnieje wiele replikowanych dość kosztownych obliczeń. Przestrzegać:

OPalgo <- function(m, p, ind1, n) {
    vcov <- matrix(0, nrow = n + 1L, ncol = n + 1)
    for (i in 0L:n) {
        for (j in i:n) {
            ## lower and upper range for the first & second multiplicand
            print(paste(c((1L + (j - i)),":",(periods - i)," 
                          ",1L,":",(periods - j)), collapse = ""))

            vcov[j + 1L, i + 1L] <- 
                sum(mat[, (1L + (j - i)):(periods - i)] *
                                mat[, 1L:(periods - j)]) /
                                    (ind * (periods - j) - 1)
        }
    }
    vcov
}

OPalgo(mat, periods, ind, n_lags)
[1] "1:70 1:70"  ## contains "1:65 1:65"
[1] "2:70 1:69"
[1] "3:70 1:68"
[1] "4:70 1:67"
[1] "5:70 1:66"
[1] "6:70 1:65"
[1] "1:69 1:69"  ## contains "1:65 1:65"
[1] "2:69 1:68"
[1] "3:69 1:67"
[1] "4:69 1:66"
[1] "5:69 1:65"
[1] "1:68 1:68"  ## contains "1:65 1:65"
[1] "2:68 1:67"
[1] "3:68 1:66"
[1] "4:68 1:65"
[1] "1:67 1:67"  ## contains "1:65 1:65"
[1] "2:67 1:66"
[1] "3:67 1:65"
[1] "1:66 1:66"  ## contains "1:65 1:65"
[1] "2:66 1:65"
[1] "1:65 1:65"

Jak widać, produkt mat[,1:65] * mat[,1:65] wykonuje się 6 razy powyżej. Jedyną różnicą między pierwszym wystąpieniem a ostatnim wystąpieniem jest to, że pierwsze wystąpienie ma dodatkowe 5 kolumn. Zamiast więc obliczać:

sum(mat[ , 1:70] * mat[ , 1:70])
sum(mat[ , 1:69] * mat[ , 1:69])
sum(mat[ , 1:68] * mat[ , 1:68])
sum(mat[ , 1:67] * mat[ , 1:67])
sum(mat[ , 1:66] * mat[ , 1:66])
sum(mat[ , 1:65] * mat[ , 1:65])

Możemy obliczyć preCalc[1] <- sum(mat[ , 1:65] * mat[ , 1:65]) jeden raz i użyj tego w pozostałych 5 obliczeniach, takich jak:

preCalc[1] + sum(mat[ , 66:70] * mat[ , 66:70])
preCalc[1] + sum(mat[ , 66:69] * mat[ , 66:69])
preCalc[1] + sum(mat[ , 66:68] * mat[ , 66:68])
preCalc[1] + sum(mat[ , 66:67] * mat[ , 66:67])
preCalc[1] + sum(mat[ , 66:66] * mat[ , 66:66])

W każdym z powyższych zmniejszyliśmy liczbę mnożeń przez 90000 * 65 = 5,850,000 i liczba dodatków wg 5,850,000 - 1 = 5,849,999 w sumie 11,699,999 operacje arytmetyczne zostały zapisane. Poniższa funkcja zapewnia to właśnie.

fasterAlgo <- function(m, p, ind1, n) {
    vcov <- matrix(0, nrow = n + 1L, ncol = n + 1)
    preCals <- vapply(1:(n + 1L), function(x) sum(m[ , x:(p - n + x - 2L)] * 
                                               m[ , 1L:(p - n - 1L)]), 42.42)
    for (i in 0L:n) {
        for (j in i:n) {
            myNum <- preCals[1L + j - i] + sum(m[, (p - n + j - i):(p - i)] * m[, (p - n):(p - j)])
            vcov[j + 1L, i + 1L] <- myNum / (ind * (p - j) - 1)
        }
    }
    vcov
}

## outputs same results
all.equal(OPalgo(mat, periods, ind, n_lags), fasterAlgo(mat, periods, ind, n_lags))
[1] TRUE

Benchmarki:

## I commented out the print statements of the OPalgo before benchmarking
library(microbenchmark)
microbenchmark(OP = OPalgo(mat, periods, ind, n_lags),
               fasterBase = fasterAlgo(mat, periods, ind, n_lags),
               RcppOrig = compute_vcov(mat, n_lags), times = 5)
Unit: milliseconds
        expr       min        lq      mean   median        uq       max neval cld
          OP 2775.6110 2780.7207 2843.6012 2784.976 2899.7621 2976.9356     5   c
  fasterBase  863.3897  863.9681  865.5576  865.593  866.7962  868.0409     5  b 
    RcppOrig  160.1040  161.8922  162.0153  162.235  162.4756  163.3697     5 a  

Jak widać, przy tej modyfikacji widzimy co najmniej 3-krotną poprawę, ale Rcpp jest nadal znacznie szybszy. Zastosujmy powyższą koncepcję w Rcpp.

// [[Rcpp::export]]
NumericMatrix compute_vcov2(const NumericMatrix& mat, int n_lags) {

    NumericMatrix vcov(n_lags + 1, n_lags + 1);
    std::vector<double> preCalcs;
    preCalcs.reserve(n_lags + 1);
    double myCov;

    int i, j, k1, k2, l;
    int n = mat.nrow();
    int m = mat.ncol();

    for (i = 0; i <= n_lags; i++) {
        myCov = 0;
        for (k1 = i, k2 = 0; k2 < (m - n_lags - 1); k1++, k2++) {
            for (l = 0; l < n; l++) {
                myCov += mat(l, k1) * mat(l, k2); 
            }
        }
        preCalcs.push_back(myCov);
    }

    for (i = 0; i <= n_lags; i++) {
        for (j = i; j <= n_lags; j++) {
            myCov = preCalcs[j - i];
            for (k1 = m - n_lags + j - i - 1, k2 = m - n_lags - 1; k2 < (m - j); k1++, k2++) {
                for (l = 0; l < n; l++) {
                    myCov += mat(l, k1) * mat(l, k2);
                }
            }
            myCov /= n * (m - j) - 1;
            vcov(i, j) = vcov(j, i) = myCov;
        }
    }

    return vcov;
}

## gives same results
all.equal(compute_vcov2(mat, n_lags), compute_vcov(mat, n_lags))
[1] TRUE

Nowe testy porównawcze:

microbenchmark(OP = OPalgo(mat, periods, ind, n_lags),
               fasterBase = fasterAlgo(mat, periods, ind, n_lags),
               RcppOrig = compute_vcov(mat, n_lags), 
               RcppModified = compute_vcov2(mat, n_lags), times = 5)
Unit: milliseconds
        expr       min         lq       mean     median         uq        max neval  cld
          OP 2785.4789 2786.67683 2811.02528 2789.37719 2809.61270 2883.98073     5    d
  fasterBase  866.5601  868.25555  888.64418  869.31796  870.92308  968.16417     5   c 
    RcppOrig  160.3467  161.37992  162.74899  161.73009  164.38653  165.90174     5  b  
RcppModified   51.1641   51.67149   52.87447   52.56067   53.06273   55.91334     5 a   

Teraz ulepszony Rcpp rozwiązanie jest około 3 razy szybsze od oryginału Rcpp rozwiązanie i około 50 razy szybciej niż pierwotny algorytm dostarczony przez OP.

Aktualizacja

Możemy zrobić jeszcze lepiej. Możemy odwrócić zakresy indeksów i / j, aby stale aktualizować preCalcs. Pozwala to na obliczenie tylko jednej nowej kolumny w każdej iteracji. To naprawdę wchodzi w grę jako n_lags wzrasta. Przestrzegać:

// [[Rcpp::export]]
NumericMatrix compute_vcov3(const NumericMatrix& mat, int n_lags) {

    NumericMatrix vcov(n_lags + 1, n_lags + 1);
    std::vector<double> preCalcs;
    preCalcs.reserve(n_lags + 1);

    int i, j, k1, k2, l;
    int n = mat.nrow();
    int m = mat.ncol();

    for (i = 0; i <= n_lags; i++) {
        preCalcs.push_back(0);
        for (k1 = i, k2 = 0; k2 < (m - n_lags); k1++, k2++) {
            for (l = 0; l < n; l++) {
                preCalcs[i] += mat(l, k1) * mat(l, k2); 
            }
        }
    }

    for (i = n_lags; i >= 0; i--) {  ## reverse range
        for (j = n_lags; j >= i; j--) {   ## reverse range
            vcov(i, j) = vcov(j, i) = preCalcs[j - i] / (n * (m - j) - 1);
            if (i > 0 && i > 0) {
                for (k1 = m - i, k2 = m - j; k2 <= (m - j); k1++, k2++) {
                    for (l = 0; l < n; l++) {
                        ## updating preCalcs vector
                        preCalcs[j - i] += mat(l, k1) * mat(l, k2);  
                    }
                }
            }
        }
    }

    return vcov;
}

all.equal(compute_vcov(mat, n_lags), compute_vcov3(mat, n_lags))
[1] TRUE

Rcpp tylko testy porównawcze:

n_lags <- 50L
microbenchmark(RcppOrig = compute_vcov(mat, n_lags),
                 RcppModified = compute_vcov2(mat, n_lags),
                 RcppExtreme = compute_vcov3(mat, n_lags), times = 5)
Unit: milliseconds
        expr       min        lq      mean    median       uq       max neval cld
    RcppOrig 7035.7920 7069.7761 7083.4961 7070.3395 7119.028 7122.5446     5   c
RcppModified 3608.8986 3645.8585 3653.0029 3654.7209 3663.716 3691.8202     5  b 
 RcppExtreme  324.8252  330.7381  332.9657  333.5919  335.168  340.5054     5 a  

Najnowsza implementacja jest teraz ponad 20 razy szybsza od oryginału Rcpp wersja i znacznie ponad 300x szybciej niż oryginalny algorytm, kiedy n-lags jest wielki.


14
2017-07-17 21:27



Miły! Dokładnie to, co miałem na myśli, kiedy powiedziałem, że można go dalej zoptymalizować. - F. Privé
Jakieś wskazówki dotyczące tego, jak włączyć lub dowiedzieć się, jak importować funkcje C ++ do swojego środowiska? próbowałem Rcpp::sourceCpp ale dostałem błędy - snoram
@snoram Czy dodałeś #include <Rcpp.h> i using namespace Rcpp; na górze pliku cpp i zastąpić # symbol komentarza z //? - Aurèle


Po prostu tłumaczenie Twojego kodu w Rcpp:

#include <Rcpp.h>
using namespace Rcpp;    

// [[Rcpp::export]]
NumericMatrix compute_vcov(const NumericMatrix& mat, int n_lags) {

  NumericMatrix vcov(n_lags + 1, n_lags + 1);
  double myCov;

  int i, j, k1, k2, l;
  int n = mat.nrow();
  int m = mat.ncol();

  for (i = 0; i <= n_lags; i++) {
    for (j = i; j <= n_lags; j++) {
      myCov = 0;
      for (k1 = j - i, k2 = 0; k2 < (m - j); k1++, k2++) {
        for (l = 0; l < n; l++) {
          myCov += mat(l, k1) * mat(l, k2); 
        }
      }
      myCov /= n * (m - j) - 1;
      vcov(i, j) = vcov(j, i) = myCov;
    }
  }

  return vcov;
}

Jest to co najmniej 10 razy szybsze niż algorytm R. Jednak czuję, że można go dalej zoptymalizować.


9
2017-07-16 20:21