在R中高效计算协方差矩阵

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

- s_baldur

---

看一下 cov()。 - Andrew Gustar

谢谢。但是如果你建议使用 cov(mat)[1:6, 1:6]，那么这有点不同...因为我不是在寻找t=1和t=2的协方差，而是在寻找一般的t和t-1之间的关系...但是如果我以不同的方式设置矩阵，也许我可以使用该函数（？）。 - s_baldur

检查?ccf函数？ - Ben Bolker

@BenBolker ccf(mat) : 仅支持单变量时间序列。 - s_baldur

如果你真的关心性能，我认为你最好的选择是使用单次遍历的Rcpp算法。我在这段R代码中唯一看到的问题是它创建了很多矩阵（副本）。 - F. Privé

Answer 1

@F. Privé的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"

如您所见，上面执行了6次产品mat[,1:65] * mat[,1:65]。第一次出现和最后一次出现之间唯一的区别是第一次出现有额外的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 解决方案比原来的 Rcpp 解决方案快大约3倍，比原始算法提供者提供的解决方案快约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  

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

- Joseph Wood

---

太好了！这正是我所想的，当时我说它可以进一步优化。 - F. Privé

有关如何将C++函数导入到您的环境中，有什么提示或学习建议吗？我尝试使用Rcpp::sourceCpp，但出现了错误。 - s_baldur

@snoram，你是否在cpp文件的顶部添加了#include <Rcpp.h>和using namespace Rcpp;，并将#注释符替换为//？ - Aurèle

Answer 2

我来翻译一下你的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;
}

这至少比R算法快10倍。 然而，我感觉它可以进一步优化。

- F. Privé