diag(X V X^T) 的紧凑/高效替代？

Question

在对线性统计模型进行预测时，我们通常有一个模型矩阵X，其中包含与我们要进行预测的点相对应的预测变量；系数向量beta；和方差-协方差矩阵V。计算预测只是X %*% beta。计算预测的方差最直接的方法是

diag(X %*% V %*% t(X))

或稍微更有效

diag(X %*% tcrossprod(V,X))

但是，这是非常低效的，因为当我们真正想要的只是对角线时，它会构造一个 n*n 矩阵。我知道我可以写一些 Rcpp-loopy 的东西来计算对角线项，但我想知道 R 中是否有一个现有的线性代数技巧可以很好地做我想要的......（如果有人想写Rcpp-loopy 对我来说作为答案我不会反对，但我更喜欢纯 R 解决方案）

FWIW predict.lm 似乎通过将 X 乘以 lm 的 QR 分解的 R 分量的倒数来做一些聪明的事情；我不确定它是否总是可用，但这可能是一个很好的起点（请参阅here）

Answer 1

按照Octave/Matlab 问题的思路，对于两个矩阵A 和B，我们可以使用AB 的nth 对角线条目将是@987654326 的乘积这一事实@ 行 A 与 nth 列 B。我们可以天真地将其扩展到三个矩阵ABC 的情况。我还没有考虑在C=A^T 的情况下如何优化，但除此之外，这段代码看起来很有希望加速：

start_time <- Sys.time()

A=matrix(1:1000000, nrow = 1000, ncol = 1000)
B=matrix(1000000:1, nrow = 1000, ncol = 1000)

# Try one of these two
res=diag(A %*% B %*% t(A)) # ~0.47s
res=rowSums(A * t(B %*% t(A))) # ~0.27s

end_time <- Sys.time()

print(end_time - start_time)

当我运行这段代码时，使用tcrossprod 似乎并没有加速结果。然而，仅仅使用 row-sum-dot-product 方法似乎已经更有效了，至少在这个愚蠢的例子中，建议（虽然我不确定）rowSums 不在返回对角线条目之前计算完整的中间矩阵，正如我所期望的那样，diag 会发生。

Answer 2

我不太确定这有多有效，

1. 找到你这样V = U %*% t(U);这是可能的，因为 V 是 cov 矩阵。
2. XU = X %*% U
3. result = apply(XU, 1, function(x) sum(x^2))

---

演示

V <- cov(iris[, -5])
X <- as.matrix(iris[1:5, -5])

使用 SVD

svd_v <- svd(V)
U <- svd_v$u %*% diag(sqrt(svd_v$d))
XU = X %*% U
apply(XU, 1, function(x) sum(x^2))
#       1        2        3        4        5 
#41.35342 39.36286 35.42369 38.25584 40.30839 

另一种方法 - 这也不会比 @davewy 的更快

U <- chol(V)
XU = (X %*% U)^2
rowSums(XU)