compact/efficient 替代 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

我不太确定这有多有效，

找到满足V = U %*% t(U)的U；这是可能的，因为 V 是 cov 矩阵。
XU = X %*% U
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)

Answer 2

沿着这个 Octave/Matlab 问题，对于两个矩阵 A 和 B，我们可以使用 nth 对角线条目 nth 的事实=14=] 将是 A 的 nth 行与 B 的 nth 列的乘积。我们可以天真地将其扩展到三个矩阵的情况，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.

compact/efficient 替代 diag(X V X^T)？

compact/efficient replacement for diag(X V X^T)?

r

linear-algebra