HC2 & HC3 optimization

[s3alfisc]

Compute $H_{ii} = diag(X \times (X'X)^{-1} \times X$) as rowSums(X * (X \times ( (X'X)^{-1})). Significantly faster for larger problems.

N <- 10000
k <- 30

beta <- rnorm(k)
X <- matrix(rnorm(N*k), N, k)
Y <- X %*% beta + rnorm(N)

tXXinv <- solve(crossprod(X))

microbenchmark::microbenchmark(
  diag(X %*% tXXinv %*% t(X)), 
  rowSums(X * (X %*% tXXinv)), 
  times = 5
)

# Unit: milliseconds
# expr         min        lq       mean
# diag(X %*% tXXinv %*% t(X)) 3725.033300 3944.9443 4319.21546
# rowSums(X * (X %*% tXXinv))    9.839701   10.3897   12.97064
# median        uq       max neval cld
# 4185.3876 4402.4162 5338.2959     5  a 
# 11.4132   16.5152   16.6954     5   b