Skuteczna obliczenie macierzy VAR KOWARIANCJA w R

Question

szukam wzrostu wydajności przy obliczaniu (auto) macierz kowariancji z pojedynczych pomiarów w czasie t z t, t-1 itp ..

W macierzy danych, przy czym każdy rząd oznacza osobę fizyczną, a każda kolumna przedstawia 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 na 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 

Answer 1

@F. Implementacja Privé to Rcpp jest dobrym miejscem startowym, 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] odbywa się 6 razy wyżej. Jedyną różnicą między pierwszym wystąpieniem a ostatnim wystąpieniem jest to, że pierwsze wystąpienie ma dodatkowe 5 kolumn. Więc zamiast computing:

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 to wykorzystać w pozostałych 5 obliczeń tak:

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, musimy zmniejszyć liczbę mnożeń przez 90000 * 65 = 5,850,000 i liczba uzupełnień o 5,850,000 - 1 = 5,849,999 dla łącznej liczby operacji arytmetycznych 11,699,999 zapisanych. 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ć z tej modyfikacji możemy zobaczyć co najmniej 3-krotny poprawę ale Rcpp jest nadal znacznie szybciej. 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 standardy:

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 wzmocnione Rcpp rozwiązaniem jest około 3 razy szybciej oryginalny Rcpp rozwiązanie i około 50x szybciej niż oryginalny algorytm dostarczonych przez PO.

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 wzrost o n_lags. Przestrzegać: tylko

// [[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 benchmarki:

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 

Ostatnio realizacja jest obecnie ponad 20x szybciej niż w wersji oryginalnej Rcpp i dobrze ponad 300x szybszy niż oryginalny algorytm gdy n-lags jest duża.

Answer 2

Wystarczy tłumaczenia 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; 
} 

To jest co najmniej 10 razy szybciej niż algorytm R. Jednak czuję, że można go dalej zoptymalizować.