0

我需要使用 Y (nxq) 和 X (nxp) 在多元线性模型中计算平方叉积矩阵的总和（实际上是该矩阵的迹线）。这样做的标准 R 代码是：

require(MASS)
require(car)

# Example data 
q <- 10
n  <- 1000
p <- 10
Y <- mvrnorm(n, mu = rep(0, q), Sigma = diag(q))
X <- as.data.frame(mvrnorm(n, mu = rnorm(p), Sigma = diag(p)))

# Fit lm
fit <- lm( Y ~ ., data = X )

# Type I sums of squares
summary(manova(fit))$SS    

# Type III sums of squares
type = 3 # could be also 2 (II)
car::Anova(fit, type = type)$SSP

这必须完成数千次，不幸的是，当预测变量的数量相对较大时，它会变得很慢。通常我只对s预测变量的一个子集感兴趣，我试图重新实现这个计算。虽然我的实现直接将线性代数转换为s= 1（下）对于小样本量（n）来说更快，

# Hat matrix (X here stands for the actual design matrix)
H <- tcrossprod(tcrossprod(X, solve(crossprod(X))), X)

# Remove predictor of interest (e.g. 2)
X.r <- X[, -2]  
H1 <- tcrossprod(tcrossprod(X.r, solve(crossprod(X.r))), X.r) 

# Compute e.g. type III sum of squares
SS <- crossprod(Y, H - H1) %*% Y

car对于大 n 仍然更快：

在此处输入图像描述

我已经尝试过Rcpp非常成功的实现，因为 R 中的这些矩阵产品已经使用了非常有效的代码。

关于如何更快地做到这一点的任何提示？

更新

阅读答案后，我尝试了这篇文章中提出的解决方案，该解决方案依赖于 QR/SVD/Cholesky 分解进行帽子矩阵计算。然而，计算所有 p = 30 个矩阵似乎car::Anova仍然比我只计算一个 (s = 1) 更快！例如 n = 5000，q = 10：

Unit: milliseconds
 expr       min        lq      mean    median        uq       max neval
   ME 1137.5692 1202.9888 1257.8979 1251.6834 1318.9282 1398.9343    10
   QR 1005.9082 1031.9911 1084.5594 1037.5659 1095.7449 1364.9508    10
  SVD 1026.8815 1065.4629 1152.6631 1087.9585 1241.4977 1446.8318    10
 Chol  969.9089 1056.3093 1115.9608 1102.1169 1210.7782 1267.1274    10
  CAR  205.1665  211.8523  218.6195  214.6761  222.0973  242.4617    10

更新 2

目前最好的解决方案是检查car::Anova 代码（即函数car:::Anova.III.mlm和后续car:::linearHypothesis.mlm）并重新实现它们以考虑预测变量的子集，而不是全部。

相关代码car如下（我跳过了检查，并简化了一点）：

B <- coef(fit)                    # Model coefficients
M <- model.matrix(fit)            # Model matrix M
V <- solve(crossprod(M))          # M'M
p <- ncol(M)                      # Number of predictors in M
I.p <- diag(p)                    # Identity (p x p)
terms <- labels(terms(fit))       # terms (add intercept)       
terms <- c("(Intercept)", terms)   
n.terms <- length(terms)
assign <- fit$assign              # assignation terms <-> p variables
  
SSP <- as.list(rep(0, n.terms))   # Initialize empty list for sums of squares cross-product matrices
names(SSP) <- terms
  
for (term in 1:n.terms){
    subs <- which(assign == term - 1)
    L <- I.p[subs, , drop = FALSE]
    SSP[[term]] <- t(L %*% B) %*% solve(L %*% V %*% t(L)) %*% (L %*% B)
}

然后只需选择术语子集。

4

1 回答 1

2

这条线和它下面的类似线H1可能会得到改进：

H <- tcrossprod(tcrossprod(X, solve(crossprod(X))), X)

一般的想法是您应该很少使用solve(Y) %*% Z，因为它与 相同solve(Y, Z)但速度较慢。我还没有完全扩展您的tcrossprod调用，以查看表达式的最佳等效公式H和H1将是什么。

您还可以查看此问题https://stats.stackexchange.com/questions/139969/speeding-up-hat-matrices-like-xxx-1x-projection-matrices-and-other-as以了解如何操作通过 QR 分解。

于 2020-10-26T10:20:50.690 回答