compact/efficient replacement for diag(X V X^T)?
Along the lines of this Octave/Matlab question, for two matrices A
and B
, we can use the use the fact that the nth
diagonal entry of AB
will be the product of the nth
row of A
with the nth
column of B
. We can naively extend that to the case of three matrices, ABC
. I have not considered how to optimize in the case where C=A^T
, but aside from that, this code looks like promising speedup:
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)
Using tcrossprod
did not appear to accelerate the results when I ran this code. However, just using the row-sum-dot-product approach appears to be a lot more efficient already, at least on this silly example, which suggests (though I'm not sure) that rowSums
is not computing the full intermediate matrices before returning the diagonal entries, as I'd expect happens with diag
.