R 通过预计算加速矩阵乘积

您的示例数据和try :

# fake data
N   = 2500
P   = 10
X   = matrix(rnorm(N*P), N, P)
W   = matrix(rnorm(P*P), P, P)
mu0 = rnorm(N)
mu1 = rnorm(N)

# precomputations
tX  = t(X)
XWX = X %*% W %*% t(X)

# functions
f_raw      = function(X, W, mu0, mu1){t(X) %*% ( ( X %*% W %*% t(X) * mu0 ) * mu1 )}
f_precomp  = function(XWX, tX, mu0, mu1){tX %*% ( ( XWX * mu0 ) * mu1 )}

更好的方法是预先计算

WX <- tcrossprod(W, X)

然后

f_better <- function (X, mu0, mu1, WX) crossprod(X, (mu0 * mu1) * X) %*% WX

笔记本电脑上的基准测试:

microbenchmark::microbenchmark(f_raw(X, W, mu0, mu1), 
                               f_precomp(XWX, tX, mu0, mu1),
                               f_better(X, mu0, mu1, WX))
#Unit: milliseconds
#                         expr        min         lq      mean     median
#        f_raw(X, W, mu0, mu1) 236.926979 238.984573 243.86319 242.212716
# f_precomp(XWX, tX, mu0, mu1) 151.190689 152.612059 156.25277 155.646093
#    f_better(X, mu0, mu1, WX)   1.113031   1.126434   1.17974   1.138207
#         uq        max neval
# 244.721163 270.477737   100
# 157.970378 182.935082   100
#   1.146352   5.130876   100

验证正确性:

ans_raw <- f_raw(X, W, mu0, mu1)
ans_precomp <- f_precomp(XWX, tX, mu0, mu1)
ans_better <- f_better(X, mu0, mu1, WX)

all.equal(ans_raw, ans_precomp)
# [1] TRUE
all.equal(ans_raw, ans_better)
# [1] TRUE

我不知道mu0和mu1在您的实际应用中是否保证为非负.如果是这样，则以下速度会快2倍:

f_fastest <- function (X, mu0, mu1, WX) crossprod(sqrt(mu0 * mu1) * X) %*% WX