R中加权最小二乘均值估计的加速逆计算答案

【问题标题】：Speed-up inverse calculation of weighted least squares mean estimate in RR中加权最小二乘均值估计的加速逆计算
【发布时间】：2014-06-13 14:15:34
【问题描述】：

我需要加快计算 R 中 WLS 中 beta 的平均估计值 - 我能够加快协方差计算 thanks to SO，现在我想知道是否还有另一个技巧可以加快平均计算（或者如果我正在做的已经足够有效）。

n = 10000
y = rnorm(n, 3, 0.4)
X = matrix(c(rnorm(n,1,2), sample(c(1,-1), n, replace = TRUE), rnorm(n,2,0.5)), nrow = n, ncol = 3)
Q = diag(rnorm(n, 1.5, 0.3))
wls.cov.matrix = crossprod(X / sqrt(diag(Q)))
Q.inv = diag(1/diag(Q))
wls.mean = wls.cov.matrix%*%t(X)%*%Q.inv%*%y
system.time(wls.cov.matrix%*%t(X)%*%Q.inv%*%y)

是否有另一个类似 wls.cov.matrix crossprod 的技巧来加速整个平均计算，还是不需要？谢谢！

【问题讨论】：

您的代码无法正常运行。此外，性能问题应包括基准和性能要求规范。
@Roland 感谢您的发现，错误已修复。在性能要求方面，我只是在寻找潜在的改进（或证明无法改进的答案），我更需要知道是否可以加快速度。

标签： r performance matrix regression

【解决方案1】：

性能的主要增益将通过沿途不使用任何 n×n 矩阵来实现。我的意思是没有 Q 矩阵，只使用它的对角线。

基于@Roland 的回答：

Qdiag = rnorm(n, 1.5, 0.3);
Q = diag(Qdiag);

wls.mean3 <- wls.cov.matrix %*% crossprod(X/Qdiag, y);
all.equal(wls.mean, wls.mean3)
microbenchmark(wls.cov.matrix %*% t(X) %*% Q.inv %*% y,
               wls.cov.matrix %*% crossprod(X, Q.inv) %*% y,
               wls.cov.matrix %*% crossprod(X/Qdiag, y),
           times=5)
#      wls.cov.matrix %*% t(X) %*% Q.inv %*% y 358050.195 363713.250 368820.818 372414.747 374824.56     5
# wls.cov.matrix %*% crossprod(X, Q.inv) %*% y  79449.856  81411.195  84616.706  85351.968  88108.62     5
#     wls.cov.matrix %*% crossprod(X/Qdiag, y)    279.092    284.867    285.252    291.796    295.26     5

【讨论】：

【解决方案2】：

在上一个问题的答案中，您被教导crossprod。再次使用该功能：

n = 1e4
set.seed(42)
y = rnorm(n, 3, 0.4)
X = matrix(c(rnorm(n,1,2), sample(c(1,-1), n, replace = TRUE), rnorm(n,2,0.5)), nrow = n, ncol = 3)
Q = diag(rnorm(n, 1.5, 0.3))
wls.cov.matrix = crossprod(X / sqrt(diag(Q)))
Q.inv = diag(1/diag(Q))
wls.mean = wls.cov.matrix%*%t(X)%*%Q.inv%*%y
wls.mean2 <- wls.cov.matrix %*% crossprod(X, Q.inv) %*% y
all.equal(wls.mean, wls.mean2)
#[1] TRUE

library(microbenchmark)
microbenchmark(wls.cov.matrix %*% t(X) %*% Q.inv %*% y,
               wls.cov.matrix %*% crossprod(X, Q.inv) %*% y,
               times=5)

#Unit: milliseconds
#                                        expr       min        lq    median       uq       max neval
#     wls.cov.matrix %*% t(X) %*% Q.inv %*% y 1019.3955 1022.1679 1022.2766 1024.540 1025.9131     5
#wls.cov.matrix %*% crossprod(X, Q.inv) %*% y  314.0622  315.3588  315.3933  317.024  317.1142     5

使用一些矩阵代数技巧可能会提高性能，但这不是我的强项。

【讨论】：