沿着Octave/Matlab问题的思路，对于两个矩阵A和B，我们可以利用AB的nth对角元素将是A的nth行与B的nth列的乘积这一事实，我们可以简单地将其扩展到三个矩阵的情况，我还没有考虑如何在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。

我不太确定这有多有效，
1.求U使得V = U %*% t(U);这是可能，因为V是覆盖矩阵。

1. XU = X %*% U
2. 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)

我最近找到了emulator::quad.diag()，它只是

colSums(crossprod(M, Conj(x)) * x)

这比@davewy的解决方案稍微好一点（尽管总体上的差异比我想象的要小）。

library(microbenchmark)
microbenchmark(full=diag(A %*% B %*% t(A)), 
               davewy=rowSums(A * t(B %*% t(A))), 
               emu = quad.diag(A,B))
Unit: milliseconds
   expr      min       lq     mean   median       uq      max neval cld
   full 32.76241 35.49665 39.51683 37.63958 41.46561 57.41370   100   c
 davewy 22.74787 25.06874 28.42179 26.97330 29.68895 45.38188   100  b 
    emu 17.68390 20.21322 23.59981 22.09324 24.80734 43.60953   100 a