r - R中var-covar矩阵的高效计算

Question

我正在寻找随着时间的推移从各个测量值计算(自动)协方差矩阵的效率增益 t与 t, t-1 ， 等等..

在数据矩阵中，每一行代表一个人，每一列代表每月的测量值(列按时间顺序排列)。类似于以下数据(尽管有更多的协方差)。

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

下面是我想出的(丑陋的)代码来获取测量/滞后测量的协方差矩阵。运行大约需要 4 秒。我确信通过移动到 data.table ，多思考而不依赖循环，我可以大大减少时间。但是由于协方差矩阵无处不在，我怀疑在 R 中已经存在一种标准(且有效)的方法来做到这一点，我应该首先了解它。

# 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。 Privé's Rcpp实现是一个很好的起点，但我们可以做得更好。您会注意到在 OP 提供的主要算法中，有许多重复的相当昂贵的计算。观察:

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"

如您所见，产品 mat[,1:65] * mat[,1:65]以上执行6次。第一次出现和最后一次出现之间的唯一区别是第一次出现有额外的 5 列。因此，而不是计算:

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])

我们可以计算 preCalc[1] <- sum(mat[ , 1:65] * mat[ , 1:65])一次并在其他 5 次计算中使用它，如下所示:

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])

在上面的每一个中，我们都减少了乘法次数 90000 * 65 = 5,850,000和增加的数量 5,850,000 - 1 = 5,849,999共 11,699,999算术运算已保存。下面的函数实现了这一点。

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

基准:

## 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  

如您所见，通过此修改，我们看到了至少 3 倍的改进，但 Rcpp仍然要快得多。让我们在 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

新基准:

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   

现在增强 Rcpp解决方案比原始解决方案快 3 倍左右 Rcpp解决方案，并且比 OP 提供的原始算法快 50 倍左右。

更新

我们还可以做得更好。我们可以反转索引 i/j 的范围以不断更新 preCalcs .这允许每次迭代最多只计算一个新列的乘积。这真的很重要，因为 n_lags增加。观察:

// [[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仅基准:

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  

最新的实现现在比原来的快 20 倍以上 Rcpp当 n-lags 时，版本比原始算法快 300 倍以上。很大。